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WSTĘP 


V drugiej połowie XIX wieku amerykański fizyk R.Wood 
odkrył ciekav/e zjawisko lokalnej kumulacji energii materiału 
wybuchowego działającego na wklęsłą metalową powłokę. W 
przeprowadzonej serii eksperymentów pokazał on, że Jeśli 
powłokę tę obciążyć symetrycznie ciśnieniem gazowych produktów 
detonacji, to z JeJ centralnych obszarów generuje się strumień 
metalu o bardzo dużej prędkości. Strumień ten może przebijać 
nawet metalowe przegrody. Energia kinetyczna cząstek strumienia 
liczona na Jednostkę masy może nawet o rząd wielkości 
przewyższać właściwą energię wewnętrzną materiału wybuchowego. 

Od tamtych czasów przyjęto nazjrwać to zjawisko, zjawiskiem 
kumulacyjnym, lub pa prostu kumulacją, choć zakres pojęciowy 
tego słowa Jest oczywiście znacznie szerszy. V pracy niniejszej 
używając terminu •‘zjawisko kumulacyjne" będę miał na uwadze 
Jego historycznie ukształtowane znaczenie, zawężone do zjawisk 
powstających przy oddziaływaniu materiałów wybuchowych na 
metalowe powłoki. V pracach 11,23 w zakres terminu "kumulacja" 
włączono również tworzenie się gazowych strumieni produktów 
detonacji w wydrążeniach materiałów wybuchowych nie pokrytych 
metalicznymi powłokami. Z punktu widzenia możliwości 
prezentowanej w pracy metody obliczeniowej, zagadnienie to 
można potraktować Jako podproblem ogólnego problemu kumulacji 
C33 . 

Koncepcje praktycznego wykorzystania efektu kumulacyjnego 
w pociskach przeciwpancernych pojawiły się w Niemczech w czasie 
I wojny światowej. Pierwszy patent na tego typu konstrukcję 
zgłoszono Już w 1914 roku. Jednakże do szerszego, praktycznego 
wykorzystania zjawiska kumulacji doszło dopiero w czasie II 
wojny światowej. Duża skuteczność broni przeciwpancernej, 
wykorzystującej efekt kumulacyjny, skłoniła wiele krajów do 
rozpoczęcia intensywnych prac eksperymentalnych 1 teoretycznych 
mających na celu zwiększenie głębokości przebicia konkretnych 
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układów i poznanie fizyki zjawiska. Pierwsza Jav/na publikacja 
na ten temat ukazała sie w 1945 roku, a jej autorami byli 
G.Birkboff, M.MacDougall, E.Pugt i G.Taylor [4]. • 

Dorobek doówiadczalno-teoretyczny prowadzonych badań, na 
który oprócz wymienionych publikacji złożyły si«ą również prace 
[5-J-12], doprowadził do powstania tzw. hydrodynamicznej teorii 
kumulacji. 

Proces generacji strumienia kumulacyjnego z weymetrznej 
warstwy ściskanej vfybuchov/o stożkowej powłoki metalowej (zwanej 
dalej wkładka kumulacyjną), opisuje s1q w ramach tej teorii 
modelem zderzających sie pod odpowiednim kątem warstw cieczy 
nieściśliwej. Parametry strumienia szacuje sią z danego 
rozkładu prędkości masoy/ej wzdłuż ramion wkładki lub v/yllczając 
go na podstawie przybliżonych ocen oddziaływania materiału 
wybuchowego na wkładką [2). Teoria ta, którą najczyściej łączy 
sly z nazwiskami Ławrientjewa i Blrkhoffa, Jest tzw. teorią 
pierwszego przybliżenia. Wyjaśniła ona podstawowy mechanizm 
tworzenia siy strumienia kumulacyjnego 1 w tych warunkach, w 
których spełnione są jej założenia, znalazła szereg 
eksperymentalnych potwierdzeń. 

Jednakże, już w początkowym okresie badań pojawiły siy 
efekty nieinożliv/e do wytłumaczenia w ramach takiego podejścia. 
Sara Lawrlentjey/ w swojej fundamentalnej pracy z 1957 roku Cl), 
podsumowując wyniki osiągnięte na bazie teorii 

hydrodynamicznej, stwierdził: "Jednakże, już teraz nagromadziła 
się dostateczna liczba faktów nie mieszczących sle w teorii i 
wymagających dla swego wyjaśnienia istotnych Jej uzupełnień". 
Zjawiska, których teoria ta wytłumaczyć nie mogła, to przede 
wszystkim; 

- istnienie krytycznego kąta rozwarcia stożkowej wkładki 
kumulacyjnej, poniżej którego nie powstaje strumień 
kumulacyj ny; 

- przemieszczanie się i zmiana parametrów strumienia 
kumulacyjnego w czasie; 


-niestacjonarne efekty towarzyszące napędzaniu wkładki i 
początkowym fazom procesu powstawania strumienia, itp. 

Dalsze prace nad teorią zjawiska kumulacji stawiają sobie 
z reguły za cel wyjaśnienie tych efektów i od połowy lat 
70-tych można pod-lelić Je na dwa dość istotnie różniące się 
kierunki. 

Do pierwszego można zaliczyć prace C13+23). Autorzy tych 
prac bazując dalej na założeniu o staćjonarności i 
nleogranlczonoścl przepływu, starają się uwzględnić dodatkowe 
elementy takie jak: ściśliwość ośrodka, lepkość, podstawowe 
charakterystyki wytrzymałościowe wkładki, itp. 

V zakresie drugiego kierunku, autorzy prac C24+36) 
wykorzystując doświadczenia tzw. fizyki komputerowej, zmierzają 
w stronę kompleksowego rozwiązania problemu, bazując na 
równaniach mechaniki ciała stałego i uwzględniając realne 
warunki początkowo-brzegowe, oraz właściwości materiałowe 
metalowych wkładek. 

Prace związane z pierwszym kierunkiem, z reguły zawierają 
analityczne formuły wygodne do stosowania w praktyce i 
umożliwiające przeprowadzanie prostych analiz wpływu niektórych 
czynników na proces kumulacji <np. ściśliwość ośrodka). 

Jednakże, daleko idące założenia upraszczające, leżące u 
podstaw tego typu rozważań, dają bardzo ograniczone możliwośćl 
wykorzystania tych wyników w praktyce. 

V drugiej grupie prac, badacze starają się z kolei 
zminimalIzować ilość upraszczających założeń i maksymalnie 
zbliżyć model matematyczno-fizyczny do realiów 
eksperymentalnych. Jednakże, informacje literaturowe odnoszące 
się do tego typu prac nie zawierają niestety, zbiorczych, 
uogólnionych wniosków, które mogłyby być pomocne przy badaniach 
konkretnych układów. Prace te zawierają z reguły przykładowe, 
ilustracyjne wyniki i wskazują w zasadzie tylko na możliwości 
określonych algorytmów, pozostających w gestii autorów. 
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Po drugie, wielość możliwych kształtów i wymiarów wkładek, 
rodzajów stosowanych materiałów wybuchowych 1 rozwiązań 
konstrukcyjno-technologlcsnych praktycznie el i minuje inożllv/ość 
sformułowania dostatecznie ogólnych 'wniosków z podawanych w 
literaturze informacji. 

V tej sytuacji, aby wykorzystać inożliv/ości Jakie daje 
modelowanie komputerowe, należy dysponować \'«rłasnym kodem 
numerycznym, który mógłby być na bieżąco wykorzystywany przy 
analizie konkretnych problemów. Praca niniejsza poświecona Jest 
właśnie opracowaniu 1 przedstawieniu możliwości takiego 
oryginalnego kodu numerycznego, umożliwiającego uzyskanie 
kompleksowych rozwiązań problemów z zakresu teorii kumulacji. 
Kompleksowość rozwiązania oznacza, iż uwzględnia ono 
nlestacjonarność 1 dwuwymiarowość zjawiska, odpowiadające 
realiom eksperymentalnym warunki początkowo-brzQgov/G i bazuje 
na półempirycznym opisie własności materiałów. Czasy obliczeń i 
obszary zajmowanej pamięci są realne w stosunku do aktualnych 
możliwości elektronicznych maszyn cyfrowych (EMO znajdujących 
się w kraju. Przykłady prezentowane w pracy zajmują obszar 
pamięci rzędu 0. 7-7-1.5 MB, a czasy ich obliczeń wahają się w 
granicach 5-^20 godzin na EMC R-60. 

V tym miejscu może zrodzić się pytanie, czy użycie tak 
skomplikav.'anego modelu fizyczno—matematyczno-numerycznego było 
nleodzovme przy modelowaniu zjawisk kumulacyjnych ? 

Okazała się, że Jeśli od wyników tego typu pracy oczekuje się 
zgodności z eksperymentem, dla wszystkich mierzalnych 
parametrów, w zakresie od kilku do kilkunastu procent, to 
użycie takiego modelu jest koniecznością. 

V pracach własnych [3, 281 pokazano, że jakościowy 
charakter zjav/lsk kumulacyjnych oddaje Już prosty, 

hydrodynamiczny model ośrodka, traktujący wkładkę kumulacyjną 
Jako tzw."ciecz sprężystą". Rezultaty znacznie lepsze 
osiągnięto zastępując ten model, modelem ciała sprężysto- 
plastycznego. Jednakże 1 w tym przypadku charakterystyki 
llQśclov/e zjawisk nie były zadowalające. Dopiero zbudowanie 
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modelu opartego na teorii sprężysto/lepko-plastycznoścl z 
uwzględnieniem efektów powstawania szczelin pozwoliło osiągnąć 
zazul er .6.ony cel, a więc zgodność wyników teoretycznych z 
eksperymentalnymi w granicach od kilku do kilkunastu procent. 
Uwagi te dotyczą przede wszystkim zjawiska, tzw. "odwrotnej 
kumulacji', które będzie omówione bardziej szczegółowo w 
rozdziale 1. V tym miejscu można Jedynie nadmienić, że zjawisko 
to polega na formowaniu pocisków z metalowych wkładek 
•stożkowych o dużych kątach rozwarcia (od 120 do 160 stopni) lub 
sferycznych o dużym promieniu krzywizny [27,374391. Zjawiskiem 
tym zainteresowano.się stosunkowo niedawno 1 Jest ono w chwili 
Obecnej intensywnie badane ze względu na swe szerokie 
możliwości zastosowań, zwłaszcza wojskowych [391. 

Teoria hydrodynamiczna nie nadaje się w zasadzie do 
11 ^c)wych analiz tego zjawiska, ponieważ w tyra przypadku nie 

zachodzi przejście metalu w fazę zbliżoną do ciekłej. Wkładka 
kumulacyjna pozostając ciałem stałym Jest tylko silnie nagrzana 
i plastycznie zdeformowana oraz zdefektowana. Na fakt znacznego 
wpływu efektów sprężysto-lepko-plastycznych na przebieg tego 
typu zjawiska zwrócono sygnalnie uwagę Już w pracach [40,411. 

Tak samo istotne znaczenie dla uzyskania odpowiedniej 
zgodności ilościowej wyników teoretycznych 1 eksperymentalnych 
ma zjawisko tworzenia się szczelin w ciałach poddanych 
dynamicznym obciążeniom. Na zjawisko to, oraz możliwości Jego 
modelowania w ramach teorii ośrodka ciągłego zwraca się 
ostatnio coraz więcej uwagi C424461. V chwili obecnej panuje 
Jednak w tej dyscyplinie duża rozbieżność poglądów i należy ją 
uznać za praktycznie całkowicie otwartą dla dalszych badań. 

Konieczność stosowania bardzo złożonego opisu 
materaatyczno-fizycznego przy próbach osiągnięcia ilościowej 
zgodności wyników teoretycznych z eksperymentalnyjiii jest Jedną 
z podstawowych przyczyn tego, że do tej pory ukazała .się 
sto-óunkowo niev/lele prac dotyczących komputerowego modeloy/ania 
zjawisk kumulacyjnych [24427, 294.31, .33-^361. 











Drugim zasadniczym powodem takiego stanu rzeczy, są znane 
trudności natury numerycznej, związane z budową odpowiedniego 
algorytmu, który dav/ałby możliwość rozwiązywania sprzężonych, 
dwuwymiarowych problemów mechaniki ciai stałych i gazów w 
warunkach skrajnie dużych deformacji. Te wielopłaszczyznowe 
wymagania powodują, że w dziedzinie tej nie można również 
wskazać na Jakiś Jeden "optymalny*' algorytm. Praktycznie każda 
z wymienionych prac ma mniej lub bardziej oryginalny charakter 
i może być traktowana również Jako próba poszukiwania takiego 
"optymalnego" algorytmu. Wydaje sią, że w chwili obecnej 
najbardziej zaawansowane wyniki osiągnięto na bazie kodu 
noszącego nazwę HEMP 1 Jego późniejszych modyfikacjach C26,2'73. 
Cechą wspólną wymienionych wyżej algorytmów Jest to, że 
wykorzystują one do budowy sieci numerycznych współrzędne 
Lagrange'a lub Eulera, a najczęściej różne modyfikacje 
mieszanych współrzędnych Lagrange'a i Eulera. Sieci te dzielą 
rozważany obiekt na komórki, w których przeprowadza się 
odpowiednie bilanse masy, pędu i energii. 

Metoda opracowana i omówiona w niniejszej pracy ma 
zasadniczo inne właściwości. Podstawową, pryncypialną cechą 
odróżniającą Ją od metod wyżej wymienionych Jest to, że nie 
twarz.y się w niej klasycznej sieci numerycznej , V metodzie tej 
operuje się Jedynie parametrami i torami wybranych na początku, 
w odpowiedni sposób punktów materialnych. Ruch i zmiany w 
czasie parametrów tych punktów wynikają z równań problemu, a 
gradienty pól w Jakich sie one poruszają wynikają z parametrów 
punktów sąsiednich, względem punktu rozważanego. Algorytm taki 
wykazuje więc większe podobieństwa do metody bicharakterystyk 
147,43] lub do tzw. "metody punktów swobodnych" [49,503, niż do 
algorytmów dotychczas stosowanych w badaniach zjawisk 
kumulacyjnych. Na wybór tego typu algorytmu wpłynęły nie tylko 
osobiste doświadczenia autora w dziedzinie komputerowego 
modelowania zjawisk z zakresu mechaniki ośrodków ciągłych 
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[43, 51-r'543, .analiza luożl Lwościj ograniczeń, zakresu 

zastosowań, stopnia złożoności, Itp.różnych metod, począwszy od 
fundamentalnych algorytmów z początku lat 60-tych [55^703. 

Ostateczną decyzję odnośnie budowy tego typu algorytmu do 
koropleksov/ego modelowania zjav/isk kumulacyjnych podjęto Jednak 
dopiero po uzyskaniu pierwszych pozytywnych rezultatów, 
przedstawionych w pracy C33. 

Obecna faza badań i uzyskane wyniki upoważniają lo 
stv/ierdzenia, że oprócz modelowania zjawisk kumulacyjnych, kod 
ten po odpowiednich adaptacjach, może być wykorzystany także do 
badania takich efekfcóv/. Jak: przebicie osłon (odpowiedni 
przykład zostanie przedstawiony w rozdziale 4>, napędzanie i 
stabilność koncentrycznych linerów, generacja fal Macha w 
układach stożkowych i walcowych C71,723, itp., a być może także 
do badania problemów wybuchowego łączenia metali. 

Na zakończenie tego rozdziału należy Jeszcze podkreślić, 
że w Polsce badania problemów kumulacji mają również 
wieloletnią tradycję i prowadzone były między innymi przez 
V.Babula 1 S.Ziębę [733, D.Smoleńskiego 1 H.Nowaka 174], 

E.Włodarczyka 1 A.Wiśniewskiego [75,763, H.Derentowicza i 
Z.Ziółkowskiego [383 oraz Z,Jopka [393. Prace te koncentrowały 
się głównie na określonych, eksperymentalnych aspektach 
problemu, a w dziedzinie teorii bazowały na wynikach teorii 
hydrodynamicznej. 

Układ niniejszej pracy Jest następujący. Po krótkim 
wstępie wprowadzającym w problematykę pracy, w rozdziale 1 
omówlćne zostaną zasadnicze elementy hydrodynamicznej teorii 
kumulacji i JeJ późniejsze uzupełnienia. Następnie, na bazie 
dostępnych informacji, krótko omówione zostaną kody numeryczne 
stosowane w badaniach zjawisk kumulacyjnych. W oparciu o te 
informacje, na zakończenie rozdziału, sformułowany zostanie cel 
niniejszej pracy. W rozdziale 2 przedstawiony zostanie układ 
równań problemu stanowiący matematyczno fizyczny model zjawska 
kumulacji. W rozdziale 3 zostaną omówione zasadnicze elementy 
kodu numerycznego służącego do rozwiązywania układu równań 
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przedstawianego w rozdziale 2. Rozdział 4 zawiera przykładowe 
rozwiązania problemów kumulacji klasycznej i “odwrotnej" oraz 
porównanie tych wyników z wynikami eksperymentu 1 teorii 
hydrodynainj,cznaJ . Na zakoiiczenie tego rozdziału zaprezentowany 
zostanie przykład modelowania przebicia osłony ciałem 
odkształcałnym, uzasadniający sformułowaną wcześniej tezą o 
znacznie szerszych możliwościach zastosowań tego kodu, nie 
ograniczonych tylko do zjawisk kumulacyjnych. Pracą zakończono 
wnioskami zamieszczonymi w rozdziale 5 oraz wykazem cytowanej 
literatury. 
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Rozdział 1 


ZJAWISKO KUMULACJI 

Rozważmy na wstąpię prosty, ideowy schemat klasycznego 
ładunku kumulacyjnego. Zbudowany Jest on z metalowej, stożkowej 
wkładki osadzonej w materiale wybucho\'fym, tak Jak to 
przedstawiono na rys.l. 



V momencie zadziałania detonatora- w materiale wybuchowym 
zaczyna slą propagować fala detonacyjna (czoło tej fali 
zaznaczono symbolicznie linią przerywaną na rys.l). Fala 
detonacyjna po określonym czasie dociera do wkładki i 
przemieszcza sią wzdłuż JeJ powierzchni. Za falą detonacyjną 
ciśnienie produktów detonacji stopniowo deformuje 1 napądza 
wkładką. W wyniku napądzanla ^Ishienty wkładki uzyskują prądkość 
o składowych radialnej i oslówej. Składowa radialna (skierowana 
w kierunku osi symetrii) powoduje gromadzenie ~lą materiału 
v/kładki wokół osi i lokalny, silny wzrost ciśnienia w tym 
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obszarze. Cl&nienle to wywału.le rozpiyw zgromadzonego przy osi 
materiału wkładki na dwie czgśct, których .masy 1 prędkości 
zależą od parametrów układu kumulacyjnego. 

Pierwsza z nich nazywa sig popularnie “zbltkiem" (po rosyjsku 
”necm”). Zawiera ona wigksza masg 1 porusza slg ze stosunkowo 
niewielka prędkością poosiową (rzgdu ułamków km/s>.T>ruga cząść. 
złożona z materiału tworzącego wewnętrzną p'Ov/iQrzchnię wkładki» 
stanowi strumień kumulacyjny, którego czoło może oslągad 
pr'^dkośó poosiową do kilkunastu km/s. Schematycznie 
przedstawiono to zjawisko na rys,2. 



fali detonacyjnej, W - prędkość napędzonych 
elementów wkładki. 

1.1. Klasyczna teoria Lawrlentj ewa-Birkhoffa oceny 
parametrów strumienia kumulacyjnego. 

Jak już wspomniano v/e wstępie, twórcami pierwszego modelu 
opisulącego zjawisko kumulacji byli niezależnie: Lawrientjew 
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[U oraz Blrkhoff, MacDougal, Pugt i Taylor C41 . Teorię tę 
nazywa się najczęściej klasyczną lub hydrodynamiczną. Zakłada 
się w niej, że wkładka kumulacyjna tworzy warstwę nieściśliwej 
cieczy idealnej napędzonej pod odpowiednim kątem do osi układu. 
Varstwa ta i“ozpływa się w otoczeniu osi na dwie części 
odpov/ladaJ ąca zbltkowi i strumieniowi ku mu lacyj nemu. Sytuację tę 
zilustrowano na rys.3. 



Rys. 3. Schemat rozpływ\j masy wkładki na strumień 1 
zbitek. 

Oznaczenia na rys.3 są następujące: 

Vs/ - prędkość ramion wkładki w układzie laboratoryjnym 
(nieruchomym) , 

^ - punkt rozdziału masy wkładki na strumień i zbitek, 
0I)» Vp* tyj - prędkość napływu masy wkładki, 

prędkość zbitka 1 strumienia w układzie 
współrzędnych związanych z punktem 
” masy elementów odpowiednio; napływających, 
zbitka i strumienia, 

OĆ ~ połówkowy kąt rozwarcia wkładki.* 
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Nie wnikając szczegółowo w problemy samego rozpływu cieczy 
t77,78], można dokonać asymptotycznego bilansu masy, pqdu i 
energii dla alemantu cieczy, który początkowo poruszał sie z 
prędkością Oi pod kątem cL do osi i miał masę a następnie 

rozpłynął sle na dwa elementy o masacb m^i oraz prędkościachi 

1 o-,-. 

Odpowiednie prawa zachowania w układzie związanym z punktem S 
Crys.3> 2 ]Qają postać? 


m t- m. 


m j m p '\y^ 


ItTol coćcl mo 




Niewiadomymi w tym układzie równań są łyi| , 1 ♦ 

Aby znaleźć czwartą, brakującą zależność, dla jednoznacznego 
określenia tych niewiadomych, można podobnie Jak w pracach [13 
1 [43 wykorzystać równanie Bernoullego, 

Dla linii prądu, która przechodzi przez zbitek, mamy; 

Pr * = c-r * 


dla linii prądu, która przechodzi przez strumień jest: 


pi + t o,- 


<1.5> i 


gdzie wielkości Pf i Pi oznaczają ciśnienia w elementach zbitka 
1 strumienia, a Cf i Cj są chwilowa dowolnymi stałymi. 
2akładając, że elementy cieczy w niezaburzonej cząścl wkładki 
są w tym samym stanie, mamy przyjąć dalej,że 

po rozdzieleniu sle cieczy na zbitek i strumień, ciśnienia w 
nich wyrównują sle, tzn. - p^* • to z równań C1.4> 1 <1.5> 
wynika, że: 
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< 1 . 6 ) 



ywzglś^ialając dalej fakt, że rozważania prowadzimy w ruchomym 
układzie współrzędnych, związanych z punktem , mamy 
dstatecźnie; 

\)rp--'xr( \ ) O'p<0 Cl.7) 

Układ równań <1,1), <1.2), <1.3> i <ł.7> pozwala wyznaczyć 
wszystkie poszukiwane wielkości: 

u-/ = 1 1 

vrp s -1 1 

fll 

% oC 

Wo co^ X 


(1.8) 

<1.9> 

< 1 . 10 ) 

Cl,li) 


Często wygodniej Jest korzystać z tych wyników 
przetransformowanych do układu laboratoryjnego Cnieruchom e go). 
jeśli w układzie tym wkładka porusza się z prędkością Vv 
normalną do swej powierzchni, to punkt rozdzielenia strumieni 
porusza się wzdłuż osi z prędkością 

-pr 

■\/ - lwi <1.12) 

~ *.no4. 

Prędkości 11 1 | W | związane są relacją; 


1 W f ^ l I oć 


<1.13) 


Stąd prędkości strumienia *V^* i zbitka w nieruchomym układzie 
współrzędnych odpowiednia wynoszą; 



I Wl 


-i -t CO40Ć 
•^'^n oi 


<1,14) 
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V'p = IVl <1.15> 

V tym miejscu należy jeszcze raz podkreślić, że 
przedstawione formuły wynikają z prostych, asymptotycznych 
bilansów masy, pądu 1 energii, zastosowanych da 
nieograniczonego, stacjonarnego przepływu nieściśliwej cieczy 
idealnej, przy dodatkowych założeniach upraszczających (stałość 
parametrów w poszczególnych cząśclach wkładkl>.Pomimo tych 
uproszczeń, teoria ta znalazła szereg potwierdzeń 
eksperymentalnych (w określonych zakresach parametrów wkładek) 

i do dziś jest stosowana, jako tzw. teoria pierwszego 
przybliżenia. Wzory <1.6), (l,9), <1,10), <1,11) pozwalają ocenić 
parametry zbitka 1 strumienia pod warunkiem, że wkładka o 
stałej grubości została napędzona do stałej prędkości pod 
stałym kątem do osi. 

V praktyce sytuacja taka najczęściej nie występuje 1 aby 
teorię tę zastosować do wkładek o zmiennych parametrach należy 
poczynić dalsze założenia upraszczające* Z reguły przyjmuje się, 

że wkładkę można podzielić na nieoddzlaływujące ze sobą 
elementy i do każdego z nich z osobna stosować wzory od <1,8> 
do <1,11) C2,163. ¥ tej sytuacji w zależności od posiadanych 
informacji lub założeń odnośnie parametrów wkładki możemy 
zbudować całą gamę różnych szczegółowych modeli. Rozpatrzmy 
bardziej szczegółowa jeden z nich, który w dalszej części pracy 
będzie można skonfrontować z rozwiązaniem numerycznym. 

Załóżmy zatem, że w chwili początkowej wkładka kumulacyjna 
o stałej grubości jest nachylona do osi pod kątem oC, a 
prędkość masowa jej elementów ma liniowy rozkład wzdłuż jej 
ramion. Sytuację tę zilustrowano na rys,4. 

Załóżmy, że rozkład prędkości j©st wzorem: 

Wiz) ~ Wo {i ~ mz.) <1.16) 

Współczynniki VVo 1 iH trzeba określić ż eksperymentu lub, tak 
jak to zostanie zrobione w dalszej części pracy, przyjąć je na 
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podstawie rozwiązań komputerowych problemu napędzania powłok 
przez materiał wybuchowy. 



Rys.4. Początkowe parametry stożkowej wkładki 

kumulacyjnej z liniowym rozkładem prędkości 
masowej 


Ha skutek zmienności prędkości 14/(z) kąt między wkładką a osią 
będzie zmieniał się w czasie ruchu wkładki. Pokazany na rys.4 
element znajdujący się w punkcie A w chwili początkowęj, w 
momencie dojścia do osi (punkt A’) jest do niej skierowany pod 
kątem , Przyjmujemy dalej, że moduł wektora IV nie ulega 
zmianie przy przejściu od A do A' oraz ,żeby dalej nie 
komplikov/ać 1 tak już mocno uproszczonego 1 przybliżonego 
podejścia załóżmy, że jest on w punkcie A' normalny do 
powierzchni wkładki « oć ^ , Z rys. 4 wynika, że; 

tq (1.17 

^ ^ A— 

ĆOSoC 


W(2:a) 















gdzie: - czaa dojścia elementu k v/łciadl:i do osi. 


dW ^ Wf: 


to Ao(. = - ^ 

T WCza-) dz. 


i po dalszych przekształceniach: 


tq Ao(.i 




Stad ostatecznie: 


5;.^ + —. 

d — ^Z.A 

«• Wo mz,A) 


^ 4- ŁOS^ 


Ponieważ punkt A jest dowolnie wybranym punktem wkładki wi^c 
współrzędną 21^ należy rozumieć jako zmienną określającą w 
chwili początkowej położenie elementów wkładki wzdłuż osi z.. 

Współrzędna Z.^ ma więc charakter współrzędnej Lagrange'a 
danego elementu. Aby określić rozkład prędkości strumienia vre 
współrzędnych Eulera, a więc wzdłuż osi w określonej chwili 
czasu, należy dokonać następującej transformacji 
przyporządkowujemy elementom Z./, ich aktualne w danej chv/ilt 
czasu położenia na osi z. , oznaczone Jako Z^(i)! 

aA(t) = 2 :a + tA • W(z-A)5mcii ■(-■V'/(z:.A)(-t-tA) (i ,24 


Oczywiście wyrażenie to ma sens dla “t po dalszych 

przekształceniach mamy; 


ZA(t)= ZA^ + tC^^et) +• WoH- WZa)’ -tA) (1.25) 

procedura wyznaczania prędkości strumienia wzdłuż osi 2. jest 
vfiQC następująca: 

- dla elementu o współrzędnej ol^reślamy ze wzoru <1«22) 
kątóJ^CzJ, a z wyrażenia (1.23) - prądkożć strumienia 

~ ze wzoru <ł.25) wyznaczamy położenie elementu strumienia 
na osi Zj w dowolnej chwili czasu Zą(^)\ 

- Jeśli tę procedurę powtórzymy dla n elementów w tej 

samej chwili czasu , to otrzymamy poszukiwany 

rozkład prędkości strumienia wzdłuż osi Z. : "VJt j. 

V ramach przyjętych założeń modelowych można wyliczyć 
jeszcze promienie elementów strumienia 1 zbitka. 

Ponieważ dla obserwatora znajdującego się w punkcie ^ zachodzi 
I a I'U'p| - ł można przyjąć* że odpowiednie masy 

elementów płynu na jednostkę długości wynoszą: 


“ TT^Tp*^ go , 

gdzie symbole 'fj-, (Tp , (T odpowiednio oznaczają: promienie 
elementów strumienia i zbitka oraz grubość wkładki w chwili 
początkowej . 

Wykorzystując wyrażenie Cl.l) dostajemy: 

'Tj*' 1- Tp^ = Ż-T-a-cT 


Z drugiej strony korzystając ze wzorów <1.10> 1 (1.11> w chwili 
gdy element znajdujący się początkowo w punkcie A osiąga punkt 
A* mamy: 


■ . i. V 

= tg 


21 











Wykorzystując wyrażenia <1.29) i <1.30>, po prostych 
przekształceniach dostajemy poszukiwany rozkład promieni 
elementów strumienia; 


Tj = /z^cTtgo4 (4- cos^)' (1.31) 

Wzór (1.31) określa Jednak tylko początkowe promienie 
strumienia w momencie Jego tworzenia sią, Istnienie gradientu 
prędkości - patrz wzór <1.23) - bidzie powodowało wydłużanie 
sle strumienia w czasie, a co za tym idzie zmniejszanie siQ 
promieni poszczególnych elementów. Jeśli przez Zą oznaczyć 
położenia początkowe elementów strumienia na osi, a przez 
aktualny promień elementu, to możemy napisać następującą 
relację wynikającą z prawa zachowania masy: 




<1.32) 


Związek między położeniem początkowym elementu na osi ^ 

Jego położeniem na wkładce 2.^Jest następujący: 

z \ » za( ”1 <1,3 

Z wyrażeń <1.32) 1 <1,33) dostajemy ostatecznie: ' 

Himo, że przedstawiona teoria zawiera wielo uproszczeń i 
modelowych założeń znalazła ona swoje eksperymentalne 
potwierdzenie w określonych zakresach parametrów wkładek. V 
następnych rozdziałach problem ten będzie szerzej rozwinięty. 


1.2. A-ktualny stan badań w ramach klasycznej teorii kumulacji. 

Szeroki zakres prac związanych z kumulacją doprowadził do 
nagromadzenia określonej liczby wyników eksperymentalnych, 
których nie można było zinterpretować w oparciu o teorię 
klasyczną, przedstawioną w głównych zarysach w poprzednim 
rozdziale. 

Wytłumaczenie tych efektów wymagało określonego 
rozbudowania tej teorii 1 określenia obszarów, w których może 
ona być stosowana. Jednym z pierwszych takich efektów, który 
nie mógł być zinterpretowany w ramach tej teorii było istnienie 
tzw. kąta krytycznego. Eksperymentalnie stwierdzono, że w 
określonych warunkach obciążenia wkładki, przy zmniejszaniu JeJ 
kąta rozwarcia o6 , poniżej pewnego kąta krytycznego strumień 
kumulacyjny nie powstaje. Tymczasem ze wzoru <1.14) wynika, że 
przy fl/—^ 0 powinno zachodzić: V^* -^oc. Wyjaśnienie tego 
zjawiska Jest następujące. Teoria hydrodynamiczna nie 
uwzględnia faktu, że w pewnych warunkach prędkość napływu 
warstwy cieczy It^jna punkt rozdzielenia ^ może być większa od 
prędkości dźwięku w materiale wkładki. 

Wprowadźmy liczbę Macha: 



gdzie C Jest lokalną prędkością dźwięku w materiale 
wkładki. 

Przy M<1 w punkcie rozdzielenia strumieni nie powstaje fala 
uderzeniowa i możemy do takiego przypadku stosować teorię 
hydrodynamiczną. Jeśli natomiast mamy przepływ z liczbą M>1 to 
w otoczeniu punktu ^ powstanie fala uderzeniowa i obraz 
zachodzących zjawisk może być Już zupełnie inny. Tak więc można 
w tym momencie sformułować pierwszy, istotny wniosek: 

klasyczna teoria kumulacji może opisywać względnie 
popravmle tylko te procesy tworzenia się strumieni, dla których 
zachodzi M<1. 
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Dla M>1 możemy mleć do czynienia z dwoma obrazami falot/ymi 
rozpatrywanego zjawiska, znan 3 'mi z mechaniki ośrodkóv/ ciągłych. 
Otóż proces odbicia fali uderzeniowej może mieć charakter 
regularny lub nieregularny w zależności od kąta zderzenia. 
Jakościowo sytuacją t^ zilustrowano na rys.5, 


a) 

M>1 

b> 


front odsuniętej 
fali uderzeniowej 

1 

7 

w 

front regularnie 
odbitej fali 
uderzeniowej 

Vo 

J 





- 

Cl 



^ s 

- ^Vi 

I 


i>i 

1 

kr 




Rys.5. Konfiguracje fal uderzeniowych w przepływie 
naddż wi Q ko wy m. 


Jak wiadomo odbicie nieregularne powstaje dla kątów > 0^ 1 . 
a regularne dla oC < Oc . W przypadku Jest oczywiste, 

że materiał v/kładkl nie może wypłynąć przez front fali 
uderzeniowej, a wi^c w tym przypadku strumieh kumulacyjny nie 
bidzie powstawał. Wartość oi ^można. wiąc określić badając kąty 

krytyczne regularnego odbicia w konkretnych materiałach, np. tak 

/ ^ 

Jak to uczyniono w pracy [79]. Kąt 0C|^^ można również określić 
na podstawie pewnych prac eksperymentalnych [10,16], 

Przykładowe wyniki wraz z innymi danymi zastaną podane w 
dalszej czyści tego rozdziału. Interesujący Jest przypadek 
kiedy M>1 i 0^> • Odsunięty front fali uderzeniowej 

(rys.5a) daje możliwość rozpłynięcia się masy wkładki 
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(po przejściu przez ten front) na struinień 1 zbitet. Okazało 
Sie. że w praktyce można takie przypadki zrealizować, ale 
strumień kumulacyjny skutkiem rozprężenia ośrodka po wyjściu z 
obszaru fali uderzeniowej ulega słabszej lub silniejszej 
fragmentacjl. Fragmenty te mogą być cząstkami ciała stałego, 
cieczy lub nawet przy oć bliskim o<^jj^mogą być parami ośrodka. W 
tym przypadku nie mamy wląc do czynienia z tworzeniem sią 
klasycznego, ciągłego strumienia kumulacyjnego. Oczywiście, 
granice miedzy litym a porozrywanym strumieniem wyznacza liczba 
«=1 i’ona też określa drugi kąt krytyczny - oć . Z (1.13) przy 
K=1 wynika,że: 



Ze wzoru <1.36> wynika bardzo praktyczny wniosek, że aby 
otrzymać Jednorodny strumień kumu 1 acyjny,prędkość Jego czoła 
nie powinna zbytnio przekraczać dwóch prędkości dźwięku w 
materiale wkładki [14,223. Zę wzorów Cl.14> 1 <1.36) mamy 

bowiem: 

■V; = c -t- ^ VV^ (1.37) 

d 

Następne ograniczenie, które należy nałożyć na teorię 
hydrodynamiczną, v/ląże się z nieuwzględnianiem przez nią 
własności wytrzymałościowych materiałów. Jest bowiem oczywiste, 
że nie powinniśmy obserwować tworzenia się strumienia, gdy 
energia kinetyczna elementów wkładki nie spowoduje plastycznego 
płynięcia materiału Cnie zostanie przekroczona granica ^ 
plastyczności Y). Wynika stąd trzeci kąt krytyczny oi = oC ^ 
który można oszacować następująco [223: 

A t ^ V Ci. 38) 

^ g W > T 

a stąd: 
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QrC COS 



(1.39> 


Rozważania te najlepie.1 bądzie podsumować posługując si^ 
przykładowymi, zbiorczymi wykresami £22], które zamieszczono na 
rys.6. 



Ciągłe strumienie kumulacyjne, które mogą być opisywane 
teorią hydrodynamiczną, powstają tylko w obszarze 1. V obszarze 
2 powstają strumienie nieciągłe (złożone z cząstek ciała 
stałego, cieczy lub pary). W obszarach 3 1 4 z omówionych wyżej 
powodów strumienie kumulacyjne nie powstają. V tym miejscu 
należy podkreślić, że parametry strumieni w obszarze 1 w miarę 
zbliżania się do kąta cx^ = ói^^będą stopniowo odbiegać od 
parametrów danych przez teorię hydrodynamiczną. Granica ta ma 
bowiem dość umowny charakter i w JeJ pobliżu własności 
wytrzymałościowe nie mogą być już pomijane. Jednakże, mimo iż 
wykresy zamieszczone na rys.6 mają nieco umowny charakter 
(krzywa ^ tG± zawierają dane dla symetrii płaskiej 
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(krzywa mają one fundamentalne znaczenie dla zrozumienia 

obserwowanych w eksperymencie zjawisk i przynajmniej 
Jakościowej ich analizy. 

Na zakoAczenle tego rozdziału należy Jeszcze zwrócić uwagę 
na inne efekty, zaobserwowane w eksperymencie, a również nie 
mieszczące się w ramach teorii hydrodynamicznej. Zauważono na 
przykład, że cylindryczne otoczki napędzone przez materiał 
v/ybuchowy do niezbyt dużych prędkości, mogą być wyraźnie 
wyhamowywane lub nawet zatrzymane przed dojściem do osi. 
Natomiast przy bardzo silnych obciążeniach zaobserwowano nawet 
parowanie wewnętrznyph ścianek otoczek skutkiem szybkiej 
zamiany energii kinetycznej w ciepło. Okazuje się, że efekty te 
mogą być wytłumaczone tylko w oparciu o złożone modele 
ośrodków, uwzględniające ich sprężysto-lepko-plastyczne 
właściwości £80,81]. Informacje te sugerują Jednoznacznie, że 
określonych rozbieżności z teorią hydrodynamiczną można 
oczekiwać w pewnych przypadkach Już na etapie zbiegania się 
wkładek do osi. Charakter tegg zbiegania może znacznie różnić 
się od przewidywanego przez model cieczy idealnej. Przykładem 
takiego trudnego przez dłuższy czas do wyjaśnienia zjawiska, 
była gfeneracja par metali w silnie obciążonych otoczkach, o 
kątach rozwarcia mniejszych niż C11,82,83] . Okazało się że 

zjawisko to wywołane Jest bardzo silnym nagrzewaniem 1 
parowaniem wewnętrznych ścianek otoczki 1 wypływem par z JeJ 
wnętrza. Jeszcze przed zajściem właściwego efektu kumulacyjnego 
(zarejestrowano prędkości par metali sięgające nawet 90 km/s). 

Przytoczone w tym rozdziale rozważania świadczą o tym, że 
w procesie kumulacji mogą występować różne zjawiska, których 
interpretacja może być trudna lub nawet niemożliwa za pomocą 
teorii hydrodynamicznej. Jedynie badania w oparciu o teorię 
sprężysto-lepko-plastyczności uwzględniającą dodatkowo efekty 
fragmentacj1, parowania, itp. są w stanie dać w miarę pełną 
interpretację eksperymentów. 
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1.3. Idea "kumulacji odvnrotnej". 

Przez dłuższy czas nie zwracano większej uwagi na zjawiska 
kumulacyjne zachodzące przy dużych kątach rozv/arcia wkładki 
kumulacyjnej. Wynikało to z faktu, że w miarą zwiększania kąta = 
zmniejsza się długość i prędkość strumienia kumulacyjnego, a co 
za tym idzie szybko spada efekty\maść Jego oddziaływania na 
przegrody (maleje głębokość przebicia). Zainteresowanie tymi 
zjawiskami zaczęło wyraźnie wzrastać w latach ?0-tych, gdy 
okazało się, że można Je wykorzystać do przebijania wszelkiego 
rodzaju osłon, ale na zupełnie innej zasadzie niż ta, na której 
oparte Jest działanie ładunku klasycznego. V literaturze 
radzieckiej [84,85] przyjęto określać to zjawisko terminem 
"kumulac.a odwrotna". Jego istotę oraz genezę tego terminu 
można wyjaśnić najprościej analizując sytuację przedstawioną na 
rys.7. 



Rys. 7.Schematyczne przedstawienie zjawisk kumulacji: 
a) - "klasycznej" ; b) - "odwrotnej". 
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lla rys.7a) przedstav/iono znany Już schemat generacji 
klasycznego strumienia kumulacyjnego. Wkładka o połówkov/ym 
kącie rozwarcia oć (w typowych ładunkach 20 -i- 30*) w wyniku 

oddziaływania na nią produktów detonacji ulega dynamicznemu 
załamaniu tworząc z osią ostry kąt oC*. V kierunku przegrody 
generuje się strumień kumulacyjny o stosunkowo nlev/ielklej 
jnasle (rzędu 20% -r 30% masy wkładki) i dużej prędkości 
czołowych elementów, rzędu 8-^10 km/s. Zbitek, zawierający 
pozostałą część masy wkładki, porusza się róvmież w kierunku 
przegrody, ale z niewielką prędkością rzędu 0,5 km/s i nie 
odgrywa praktycznie żadnej roli w procesie jej przebijania. 

Natomiast na rys.7b) przedstawiono Ideowo proces 
oddziaływania produktów detonacji na wkładkę o dużym (o^ =: 60* 4- 
80*), początkowym kącie rozwarcia. W wyniku tego oddzlałyv/anla 
dynamiczny kąt załamania v/kładkl oć* może być większy od kąta 
prostego. Według metodologii klasycznej teorii kumulacji (wzory 
(1.8) 4 - (I.ID) będziemy mieli w tym przypadku następującą 
sytuację: v/ kierunku przegrody będzie przemieszczać się 
najpierw stosunkowo duża część masy wkładki (=? 70% 4- 80%) z 
dość znaczną prędkością rzędu 2 4-3 km/s, a za nią podążał 
będzie niewielki strumień o kilka razy mniejsi 2 j prędkości. 
Koźna więc powiedzieć, że w stosunku do kumulacji klasycznej 
"odvnrócone" zostały role zbitka i strumienia. Elementem 
szybszym i decydującym o przebiciu jest teraz zbitek, a nie 
strumień. Zgodnie z taką umowną klasyfikacją o podziale zjawisk 
kumulacyjnych na "klasyczne" i "odvfrotne" decyduje dynamiczny 
kąt załamania wkładki . Granicę rozdziału wyznacza kąt cC - 
90*. 

Okazało się Jednakże, że taka klasyfikacja i nazewnictwo 
nie są adekwatne do zjawisk obserwowanych w eksperymentach. 
Wynika to stąd, że przy dużych kątach cL decydujący wpływ na 
przebieg zjawiska mają efekty wytrzymałościowe. Model 
klasyczny, w którym wkładkę traktuje się Jak warstwę cieczy 
idealnej, nie może w takiej sytuacji oddać ilościowych ani 
nawet Jakościowych charakterystyk procesu. Efekty 
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wytrzyioaiościowe powodują bowiem szybkie wyróvmywaaie, 
stosunkowo niewielkich w tym przypadku gradientów prędkości 
poosiowej i uniemożliwiają rozdzielenie przepływu na strumień i 
zbitek- Praktycznie cała masa zdeformowanej wkładki, tworząc 
Jedno ciało przemieszcza slq w kierunku przegrody. Kształt tego 
ciała <poclsku) zależy od konkretnych charakterystyk całego 
układu kumulacyjnego. Z tego powodu coraz czyściej proponuje 
sie zastąpienie terminu '‘kumulacja odwrotna” sformułowaniem 
"wybuchowe kształtowani© pocisków", 

V chwili obecnej ładunki do wybuchowego kształtowania 
pocisków znalazły Już szereg praktycznych zastosowań, mimo iż 
głębokość przebicia pancerza takim pociskiem Jest od 5 do 10 
razy mniejsza od klasycznego ładunku kumulacyjnego o tej samej 
średnicy wkładki. Ten niedostatek jest jednak rekompensowany 
dużą masą pocisku i prawie o rząd większą średnicą krateru jaki 
tworzy on w pancerzu. Jest więc oczywiste, że w tych 
przypadkach, w których nie Jest wymagana zbyt duża głębokość 
przebicia, efekt niszczący takiego pocisku będzie o wiele 
silniejszy niż klasycznego strumienia kumulacyjnego, 

Z drugiej strony wiadomo, że klasyczne ładunki kumulacyjne 
muszą być detonowane w dość precyzyjnie dobranej odległości od 
pancerza, aby osiągnąć maksymalną głębokość przebicia. 
Zdetonowanie takiego ładunku na jakiejkolwiek przeszkodzie 
przed właściwym pancerzem powoduje, że strumień zbyt wcześnie 
fragmentuje się i rozprasza radialnie, Prowadzi to do szybkiego 
spadku głębokości przebicia. Natomiast wybuchowo ukształtowany 
pocisk o dużej masie i wyrównanej prędkości poosiowej może 
Jednakowo skutecznie razić cele znajdujące się w różnych 
odległościach od miejsca detonacji. Odległości te mogą wynosić 
dziesiątki, a nawet setki metrów. Na torze swego lotu pocisk 
taki maże przebijać nie tylko pojedyńcze, ale i kaskadowo 
ustawione przegrody. Te właściwości wybuchowo kształtowanych 
pocisków zdecydowały o ich licznych zastosowaniach np. w minach 
dennych, burtowych 1 Innych konstrukcjach, zwiększając znacznie 
ich kierunkową siłę i zasięg rażenia. 
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Pierwszą minę przeciwpancerną wykorzystującą ten efekt 
zaprezentowano w RFN w 1962 r. ¥ latach 70-tych rozpoczęto 
bardzo Intensywne badania tych zjawisk, które w latach 80-tych 
weszły na Jakościowo wyższy etap, wspomagane złożonymi modelami 

V Polsce badania eksperymentalne zjawiska "odwrotnej 
kumulacji'* zostały zainicjowane w WAT i prowadzone we 
współpracy z IFPiLM (patrz np. prace [38,39]>. Wyniki tych 
badań wpłynęły inspirująco na podjęcie tematu niniejszej pracy, 

1,4. Informacje literaturowe dotyczące modelowania 
komputerowego zjawisk kumulacyjnych. 

Silne założenia' upraszczając© leżące u podstaw klasycznej 
teorii, kumulacji oraz praktycznie brak możliwości ich Istotnego 
złagodzenia przy, próbach znalezienia rozwiązań analitycznych, 
zrodziły potrzebę poszukiwania Innych sposobów rozwiązania tych 
problemów, 

% pomocą przyszły tu metody fizyki komputerowej, intensywnie 
rozwijane zwłaszcza w USA i ZSRR od końca lat 50-tych. 

Jednakże, mimo zgromadzenia dużej liczby prac z teJ dziedziny, 
zagadnienia kumulacji dość długo nie mogły doczekać się choćby 
uproszczonych rozwiązań. Modelowanie komputerowe zjawisk 
kumulacyjnych wymaga bowiem rozwiązania szeregu bardzo trudnych 
problemów, o których wspomniano już we wstępie, a które 
dokładniej zostaną przedstawione w punkcie 1,6, 

Pierwsza praca t24] dotycząca modelowania komputerowego 
zjawisk kumulacyjnych ukazała się dopiero w 1976 1 w dalszym 
ćią.gu informacji o takich modelach jest niewiele. W chwili 
obecnej informacje odnośnie modelowania komputerowego 
•kumulacji klasycznej", można znaleźć w pracach [24-ł-263 , a 
•kumulacji odwrotnej" w C27,29-^36] , 

Prace dotycząc© "kumulacji odwrotnej" zaczęły ukazywać się 
praktycznie dopiero w połowie lat 80-tych. V ostatnim czasie 
Ukazało się również kilka pozycji [86,873 w specjalistycznych 
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czasopismach wojskowych, z których wynika, że rozpoczęto Już 
stosowanie niektóry- h isodeli do projektowania konkretnych 
układów kumulacyjnych, W oparciu o model, który bidzie 
zaprezentowany w pracy niniejszej, badania wspierające 
projektowanie układów kumulącyJnych rozpoczęto w 1987 r. 

Na bazie dostępnych informacji literaturowych można w 
chwili obecnej zbudować pewną ogólną klasyfikacje metod 
komputerowych, które zostały wykorzystane w badaniach zjawisk 
kumulacyjnych, W zależności od typów wykorzystywanych schematów 
numerycznych, są to: 

- metoda współrzędnych Lagrange'a [27,33]; 

- metoda współrzędnych Eulera [24,25,26,31]; 

" metody mieszanych współrzędnych Eulera i Lagrange*a 
[29,30,32,34,35]; 

^ metoda elementów skończonych [36], 

Przed bardziej szczegółowym przeanalizowaniem najistotniejszych 
wyników tych prac można ogólnie następująco scharakteryzować te 
grupy algorytmów: 

i. Metoda współrzędnych Lagrange’a, 

V algorytmach tego typu numeryczna sieć porusza sle i 
deformuje wraz z ośrodkiem. Większość prac z tej dziedziny 
bazuje na fundamentalnym algorytmie M. Wilkinsa [55],, Podstawową 
zaletą takiego algorytmu Jest to, że umożliwia on wygodne 1 
dokładne modelowanie ruchomych granic rozdziału ośrodków i 
linii brzegowych. Po drugie, bez, względu na zachodzące 
deformacje, zachowuje się w trakcie rozwiązywania problemu 
stała liczba węzłów” sieci numerycznej. Po trzecie, sieć 
poruszająca się wraz z ośrodkiem pozwala na uniknięcie szeregu 
problemów związanych z transportem strumieni masy, pędu czy 
energii przez granice siatek. Podejście takie ma Jednak Jeden 
bardzo poważny mankament, Otóż, Jeśli w wyniku ruchu ośrodka 
sieć Lagrange’a zaczyna się silnie deformować, wówczas 
pierwotnie regularne komórki sieci zaczynają przybierać 
stopniowo coraz bardziej nieregularne kształty. Prowadzi to z 
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cjjasem do utraty dokładności obliczeń, a utrzymanie stabilności 
wymaS^ szybkiego zmniejszania kroku czasowego, ¥ takim momencie 
ńal«^y obliczenia przerwać, albo - Jeśli zależy nam na analizie 
dalszych faz zjawiska - dokonać rekonstrukcji siatki [883, 

24 Xetoda współrzędnych Eulera. 

¥ przypadku bardzo silnych deformacji ośrodka metoda 
Lagrange'a może okazać się zbyt niewygodną w stosowaniu. Lepsze 
rezultaty można wówczas osiągnąć stosując algorytmy 
wykorzystujące współrzędne Eulera. V algorytaaach tych stosuje 
się nieruchomą, regularną sieć współrzędnych, przez którą 
przepływa rozważany ośrodek. V ten sposób można dość prosto 
uniknąć komplikacji związanych z dużymi deformacjami ośrodka, 
5 iXe w zamian za to pojawiają się inne, dość poważne problemy. 
Granice ośrodka przemieszczające się na tle nieruchomej siatki 
tworzą zmieniające się w czasie, nieregularne komórki, w 
których obliczenia należy prowadzić za pomocą specjalnych 
algorytmów wykorzystujących elementy metody Lagrange'a 
(wprowadza się np. tzw. ”cząstki znaczone** poruszające się vnraz 
z ośrodkiem C61]>, Znaczne trudności wynikają również przy 
obliczaniu transportu strumieni przez nieruchome oczka statki, 

Z reguły proces obliczeniowy rozbija się na kilka etapów, np, w 
pierwszym etapie nie uwzględnia się członów konwekcyjnych i 
wylicza tylko przybliżone wartości wynikające z pracy sił 
ciśnienia, a w drugim etapie korzystając z tych </artośoi 
wylicza się strumienie między poszczególnymi komórkami siatki 
[69,70], Zaawansowane i dokładne algorytmy [69] zawierają z 
reguły więcej rozbudowanych etapów obliczeń i korekcji, i 
prowadzą do bardzo złożonych struktur obliczeniowych. 

3, Metoda mieszanych współrzędnych Eulera i Lagrange*a. 

V wielu zagadnieniach, ze względu na wyżej wymienione wady 
obu metod, proponuje się używanie algorytmów mieszanych, 
wykorzystujących ich elementy. 
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z tej grupy metod w zagadnieniach kumulacji znaiasrły do taj 
pory zastosowanie głóymie; 

^ metoda "cząstek w sieci" (PIC); 

- metoda "dużych cząstek"; 

~ metoda "ruchomych siatek". 

V metodzie "cząstek w sieci" (PIO kontinualny ośrodek 
^astepuje slą zbiorem cząstek, które unoszą ze sobą masą, pąd i 
energie ośrodka 1 przemieszczają sle na tle nieruchomej, 
eułerawsfciej siatki współrzędnych. Zbiór poruszających sle 
cząstek tworzy analog drugiej, lagrange-owsktej sieci. Stec 
nieruchoma s^ży do wyznaczania gradientów pól, w których 
poruszają sie cząstki ośrodka. Niewątpliwą zaletą tej metody 
Jest możliwość modelowania złożonych ruchów wieloskładnikowych 
ośrodków [53,61,1301. Wadami tej metody są natomiast wysokie 
wymagania odnośnie mocy obliczeniowej EMC (pamięć 1 szybkość 
obliczeń) oraz trudne do wyeliminowania fluktuacje występujące 
Wówczas gdy w jakiejś komórce znajdzie sią zbyt mała liczba 
dyskretnych cząstek. Oczywi.ścle, można starać sie przezwyciężyć 
te trudności różnymi sposobami [891 <np. przyjmując określone 
"rozmycie dyskretnych cząstek") ale i tak, wyniki tej metody 
uważa sle za lepiej oddające Jakościowe niż ilościowe 
Charakterystyki zjawisk. 

Ta podstawowa wada metody PIC była punktem wyjścia do 
budowy tzw. "metody dużych cząstek". Nie postuluje sle w niej 
rozbicia ośrodka na dyskretne cząstki, ale traktuje masą 
zawartą w komórce Eulera Jako Jedną "dużą cząstką" 1 rozważa 
Sie ciągłe strumienie tych "dużych cząstek" przez granice 
siatki [671. V tym miejscu nasuwa sle uwaga, że metoda ta 
wycholząc z lagrangeowsko-eulerowsklej metody PIC staje sie 
praktycznie bardzo bliska metodom określanym przez Innych 
autorów Jako eulerowskle (por. prace C671 i 1701). 

Metoda "ruchomych -latek" Jest również próbą połączenia i 
wykorzystania zalet mc Lagrange-a 1 Eulera. Wprowadza sle w 
niej ruchoma, krzywoliniowe współrzędna, które wygodnie Jest 
wybierać tak. aby linie brzegowe oddawały kształt granic 
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rozważanego cia.la. Z reguły Jedna ze współrzędnych ma charakter 
^ ^^półrzędnej Eulera (nieruchoma)» a druga ma charakter 

rządnej Lagrange'a (przemieszcza sią wraz z ośrodkiem) 
130,60]. Konkretny wyhór sieci Jest zdeterminowany charakterem 
ruchu ośrodka. Optymalny dobór takiej sieci jest Jednak w 
praktyce bardzo skomplikowany £90], a złożoność tego 
zagadnienia Jest porównywalna ze złożonością całego wyjściowego 
problemu, 

V ostatnich latach pojawił stą również przykład 
zastosowania metody elementów skoilczonych [361 do modelowania 
zjawiska "odwrotnej kumlacji". Wychodząc z zasad mechaniki w 
sformułowaniu wariacyjnym metoda ta dość istotnie różni sią od 
typowych metod różnic skoilczonych. Operuje sią tu bowiem 
ciągłym przyblźeniem rozwiązania ścisłego i dąży do tego aby 
było ono Jak najlepsze w ramach przyjętego systemu funkcji 
bazowych, Osiąga sią to poprzez minimalizacją odpowiednich 
funkcjonałów, dochodząc w efekcie do systemów odpowiednich 
równań różnicowych. Zaletą takiego podejścia Jest duża swoboda 
vf sposobie dyskretyzacj 1 ośrodka i co za tym idzie, ułatwione 
modelowanie układów o bardzo złożonych geometriach. Najlepsze 
rezultaty metoda ta daje dla równań typu eliptycznego. Dla 
równań hiperbolicznych JeJ efektywność staje sią dyskusyjna. 
Przede wszystkim prowadzi ona w ścisłym sformułowaniu do 
niejawnych schematów różnicowych, które są kłopotliwe przy 
rozwiązywaniu i wprowadzają niepożądaną dyfuzją numeryczną 
prowadzącą do niefizycznego "rozmywania" frontów falowych. 
Oprócz tego, w zagadnieniach z ruchomymi brzegami Jest sens 
stosować metodą do równań zapisanych we współrządnych. 

Lagrange*a. Rozwiązanie takie bądzie wiąc obarczone tymi samymi 
ograniczeniami, które są charakterystyczne dla wyżej omówionych 
. .. metod Lagrange' owskich. 

Zakres, objątość i cel niniejszej pracy nie pozwalają na 
ogólne scharakteryzowanie nawet cząści algorytmów, ich różnych 
wersji i konkretnych sposobów realizacji, poszczególnych zadań, 
które dotyczą ogólnie rozumianej mechaniki ośrodków ciągłych. 

















Nawiązano Jedynie bardzo ogólnie do tych grup metod, dla 
których znaleziono p-zykłady ich zastosowań w .badaniach zjawisk 
kumulacyjnych. Aktualność problematyki i wzrastająca szybko 
ilość publikacji na ten temat wskazują, że należy oczekiwać 
zarówno doskonalenia istniejących metod v/ tej dziedzinie jak 
również pojawienia się nowych, oryginalnych algorytmów. 

Po ogólnym scharakteryzowaniu metod obliczeniowych, na 
zakończenie tego rozdziału, wydaje sią celowym krótkie 
omówienie podstawowych rezultatów prac dotyczących modelowania 
komputerowego zjawisk kumulacyjnych. Jak Już wspomniano, za 
pierwszą z tych prac można uznać prac^ Chou, Corleone i Karppa 
£24]. V pracy tej przyjęto, że płaska lub stożkowa wkładka 
kumulacyjna została natychmiastowo napędzona do określonej 
prędkości (normalnej do powierzchni wkładki) i w chwili 
początkowej Jest nachylana do osi pod stałym kątem. Wprawdzie 
takie warunki nigdy nie występują w eksperymencie, ale celem 
tej pracy była nie analiza realnych wkładek kumulacyjnych, a 
Jedynie weryfikacja kryteriów powstawania strumienia 
kumulacyjnego. Przyjmowano wląc warunki początkowe tak aby 
znaleźć sie w obszarach 1,2 lub 3 przedstawionych na rys.6. 
Wykorzystali przy tym realne, półemplryczne charakterystyki 
materiałowe. Praca potwierdziła poprawność ocen, z których 
korzysta s1q przy formułowaniu kryteriów powstawania ciągłych 
lub pofragmentcwanych strumieni kumulacyjnych. Szersze badania 
tego typu mogą być wykorzystane do budowy wykresów 
przedstawionych na rys.6 dla szerokiej klasy materiałów, bez 
wprowadzania dodatkowych założeń upraszczających. 

Z kolei w pracy £253 przedstawiono ide^ konstrukcji 
eulerowsklej wersji kodu numerycznego, za pomocą którego 
uzyskano rozwiązanie dla wstępnej fazy procesu tworzenia si^ 
strumienia. W rozwiązaniu tym uwzględniono Jednocześnie 
detonację materiału wybuchowego oraz deformację i napędzanie 
wkładki. Jednakże zaprezentowanie Jednego wyniku dla bardzo 
wczesnej fazy zjawiska nie pozwala ocenić czy kod ten daje 


36 


vrynlkl zgodne z eksperymentem, ani czy pozwala na kontynuowanie 
obliczeń dla dalszych faz zjawiska. 

Praca £263 stanowi próbę zastosowania kodu typu Lagrange’a 
do modelowania zjawiska kumulacji klasycznej. Jak należało 
oczekiwać kod tego typu ze względu na swe ograniczenia może być 
użyty Jedynie do analizy wstępnych faz procesu, a praktycznie 
Jedynie do uzyskania informacji o stanie wkładki (kształt 1 
prędkość) po JeJ napędzeniu przez produkty detonacji. Autorzy 
pracy w oparciu o tę informację analizują proces tworzenia się 
i propagacji strumienia w/g metodologii klasycznej teorii 
kumulacji analogicznej do przed.stawion.eJ w rozdziale 1.1. 
Uzyskano przy tym dość dobrą zgodność niektórych parametrów 
strumienia z danymi eksperymentalnymi, pomimo przyjęcia 
arbitralnych założeń zv/iązanych z wykorzystaniem wzorów 
podanych w rozdziale 1. Praca ta świadczy dobitnie o wielkich 
trudnościach Jakie muszą być pokonane, aby w ramach Jednego 
rozwiązania uwzględnić zarówno proces napędzania wkładki Jak 1 
proces tworzenia się strumienia kumulacyjnego. Przykładów 
takich kompleksowych rozwiązań skonfrontowanych z eksperymentem 
do chwili obecnej nie znaleziono w literaturze. 

Zapewne, z uwagi na znacznie mniejsze deformacje oraz 
korzystniejsze (z punktu widzenia możliwości analiz 
numerycznych) relacje wymiarów wkładek, liczba prac 
poświęcanych zagadnieniu "kumulacji odwrotnej" Jest Już w 
chwili obecnej znacznie większa £27,29-5-363. 

Jednakże, wyniki przedstawione w tych pracach, w Jeszcze 
większym stopniu niż dotyczyło to kumulacji klasycznej, noszą 
charakter Jedynie przykładowych ilustracji możliwości metod 
obliczeniowych. Nie mają one natomiast charakteru analizy 
zjawiska w zależności od zmian określonych parametrów. 

Po drugie, z reguły nie podaje się większości parametrów 
rozważanych układów i po trzecie, najczęściej v/yniki te dotyczą 
początkowych faz procesu bez określenia ostatecznego kształtu 
uformowanych pocisków. Na osobną uwagę zasługują w tej’ 
dziedzinie prace C27,33,86,873. 
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2 prac tycłi nie można również wyciągnąć wniosków odnośnie % 
zakresu i rezultatów prowadzonych analiz* ale pozwalają one J 
przynajmniej na stwierdzenie, że w przypadkach niezbyt dużych § 
def ormacJ i"® uzyskano dobrą zgodność cyników teoretycznych z m 

eksperymentalnymi. Zgodność ta dotyczy również późnych faz 1 
zjawiska, w których można już mówić o uformowaniu sią pocisku* M 
Ka zakończenie tego rozdziału, należy odnotować jeszcze, żeM 
w zakresie niezbyt dużych deformacji oraz dla kształtów wkładek! 
nie wymagających zbyt dużej liczby podziałów przedstawiono już M 
pierwsze rozwiązanie problemu sformułowanego w trzech M 

niezależnych zmiennych przestrzennych [363. Badania tego typu || 
mogą być pomocne w ocenach wpłyvAj eksperymentalnych asymetrii M 
na proces tworzenia si^ wybucho\ifo kształtowanych pocisków. | 

1.5. Cel pracy. J 

Z analizy przedstćiwionego laaterlału wynika, że klasyczna J 
teoria kumulacji nie opisuje adekwatnie do rzeczywistości J 

kompleksu zjawisk towarzyszących napędzaniu 1 deformacji powłok! 
<wkładek) metalowych o złożonej geometrii w wyniku J 

oddziaływanla na nie produktów detonacji. Kie uwzględnia ona li 
fIzyko-mechanicznych właściwości wkładek, jak również procesów 
związanych z detonacją 1 właściwościami materiałów wybuchowych. H 
Dane literaturowe z zakresu fizyki komputerowej, które 
uwzględniają wymienione właściwości wkładek i materiałów {[ 

wybuchowych, podawane są w formie sygnalnej i nie można Ich 0 

wykorzystać do analizy konkretnych układów kumulacyjnych, 

Wynikła zatem potrzeba opracowania nowego mateiaatycznO’“ U 
fizycznego modelu układów kumulacyjnych wraz z odpowiednim ^ 

kodem numerycznym, umożliwiającym komputerowe symulowanie li 

dynamiki procesów napędzania i programowanego kształtowania 
metalowych powłok, w pełnej zgodoścl z warunkami eksperymentu, 'i\ 
Dla spełnienia wymienionych wymagań model taki powinien i\ 
zawierać? I? 

ogólne równania gazodynamiki i półempiryczne równania 


^tanu dla opisu przestrzennych procesów detonacji 
loateriałów wybuchowych; 

^ równania teorii sprężysto/lepko-plastyczności wraz z 
półempirycznymi charakterystykami materiałowymi dla 
opisu procesów deformowania i napędzania metalowych 
wkładek kumulacyjnych. 

Algorytm obliczeniowy powinien natomiast umożliwiać; 
prowadzenie obliczeń w warunkach skrajnie dużych 
odkształceń i niejednorodności zdeformowanej wkładki, 
łącznie z możliwością JeJ fragmentacji; 

- stawianie warunków brzegaV7ych na ruchomych 1 krzywoliniowych 
powierzchniach; 

- zszywanie rozwiązań na kontakcie metal - produkty detonacji; 

- prowadzenie analiz efektów kumulacyjnych dla wkładek o 
różnych kształtach i wymiarach; 

- otrzymanie stabilnych rozwiązań dla dostatecznie długich 
czasów (w badanym problemie do czasów oderwania się 
strumienia kumulacyjnego od zbitka lub ostatecznego 
uformowania się pocisku w procesie '^kumulacji odwrotnej”). 

Skonstruowanie i przedstawienie takiego kompleksowego 

modelu wi-az z przykładami Jego zastosowań skonfrontowanymi z 

eksperymentem jest celem niniejszej ej pracy. 
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Rozdział 2 


]UTEKATYCSSO~ SFORMUŁOWANIE PROBLBKTJ 

V rozdziale tym wictnie zastaną równania, którymi 

opisano procesy detonacji ‘Skrmdensowanych .materiałów 
wybuchowych <KV> oraz deformacji i napędzania metalicznych 
wkładek kumulacyjnych. Róvmania te opisują zjawiska 
niestacjonarne i dwuwymiarowe przestrzennie. W zagadnieniach 
kumulacji najwygodniej Jest stosować równania zapisane w 
symatrii Cylindrycznej < r , ^ * O ) 1 w takiej też 

postaci bqdą one używane w niniejszej pracy. Można tylko dodać 
że metoda numeryczna, która bidzie omówiona w nastąpnym 
rozdziale, pozwala na bardzo proste i szybkie przejście do 
rozwiązywania zagadnień w symetrii płaskiej ^ ^ ^ ^ “ 0 ), 

Podstawowe równania gazodynamiki i mechaniki ciał 
odkształcałnych nie bądą w pracy wyprowadzane, ani szeroko 
komentowane, gdyż są one znane z bardzo wielu publikacji. 
Bardziej szczegółowym komentarzem zostaną natomiast opatrzone 
te elementy matematyczno-fIzycznego opisu, które są wynikiem 
najnowszych opracowań lub mają charakter własnych badań. 

Dotyczy to głównie efektóv/ lepkich oraz zjawisk tworzenia się 
szczelin w ciałach stałych. 

Przedstawione w poprzednim rozdziale problemy Jakie 
występują przy analizie zjawisk kumulacyjnych oraz dążenie do 
uzyskania maksymalnej zgodności między wynikami teoretycznymi 
eksperymentalnymi spowodowało, że model matematyczno-fizyczny 
był w trakcie pracy ciągle rozbudowywany 1 uzupełniany. V swej 
aktualnej postaci Jest on modelem dość ogólnym ł uwzględnia 
najnowsze wyniki badań dotyczące własności materiałów w 
warunkach dynamicznych obciążeń. Zaletą takiego ogólnego opisu 
Jest to, że pozwala on w ramach jednego modelu numerycznego na 
uzyskanie zgodnych z eksperymentem rozwiązań szeregu złożonych 
zagadnień początkowo-brzegowych (kumulacja klasyczna, kumulacj 
odwrotna, przebijanie przegród pociskami itp.>. 
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2.1' Opis procesu detonacji materiałów viybuchowvch. 


Ko spisu tprocesu detonacji będziemy stasować teorię 
hydrodynaiuiczną. Zgodnie z tą teorią zachowanie się produktów 
detonacji można opisać znanym układem równań: 

4t ^ o < 2.1 


t> 

S ti^ 


przy następujących oznaczeniach: 


fC * X -współrzędne przestrzenne (symetria cylindryczna) 
g - gęstość, 
p - ciśnienie, 

AL I ^ ~ składowe wektora prędkości ruchu produktów wybuchu tO' 
odpowiednio wzdłuż osi IT i Z , 

^ - energia wewnętrzna na Jednostkę masy. 


Pochodna substancjalna -gpę- oraz wyrażenie div UT , są 
zdefiniowane następująco: 


i? * Ir » 7 


Układ równań (2,1) do (2,4) uzupełnia równanie stanu 
produktów detonacji, wiążące ciśnienie z gęstością i energią 
wewnętrzną. 

V literaturze można znaleźć bardzo wiele różnych postaci 
równań stanu dla szerokiej gamy materiałów wybuchowych. 
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Syntatyczny przegląd tych rówia* z os ta 3: podany, 
C9i] . 



Przy)£ładawe 2 'azwiązania zamieszczone w pracy niniejszej | 

bazowały początkowo na równaniu stanu termoplastycznego K¥ o | 
postaci; 

p= 03 » + cTsł (2,?)i 

I 

gdzie; C I ^ " stałe empiryczne. 

V pó^żniejszych pracach wprowadzono Jednolitą i bardzo dokładniel 


opracowaną postać równania stanu przez tzw* zespół EOS z 
laboratorium Los Alaraos, Równanie to oznaczono symbolem JVL 
(Jones, Wllkins, Lee^ i podano wartości odpowiednich 
współczynników dla szerokiej gamy różnych materiałów 
wybuchowych. 

Równanie stanu JVL produktów detonacji ma następującą 
ogólną postać; 



Wartości stałych współczynników A , 5 . , R. 2 , i 

typowych MV zebrano w tabeli 1. I 

W rozważanych zagadnieniach kumulacji, wymiary układów są | 
o kilka rzędów wielkości wląksze od wymiarów | 

charakterystycznych strefy reakcji chemicznych w fali 
detonacyjnej. Pozwala to na stosowanie hydrodynamicznej teorii | 
detonacji i modelowanie frontu fali detonacyjneJ powierzchnią 
silnej nieciągłości. j 

V Wielu pracach [20,32,44,92] sprawdzono, że daje ono dobre i 
rezultaty i co również Jest bardzo ważne - nie wymaga zbyt i 

geotych sieci numerycznych. Polega ono na tym, że zakłada sle 
iż front fali detonacyjnej jest ruchomą powierzchnią silnej | 
nieciągłości o zadanym kształcie i parametrach odpowiadających | 
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Tabela i. Współczynniki dla równania stanu typu JVL. 


Uazwa 

UW 

S- 

i g/cra"^] 

A 

tKbar] 

3 

[ Kbar} 



cT 

komp.A"3 

1.650 

6.. 113 

0,1065 

4.40 

l. 20 

0. 32 

Komp. B 

1.717 

6.242 

0.0768 

4.20 

1. 10 

0.34 

Komp. C-4 

1.601 

6. 098 

0.1295 

4. 50 

1.40 

0. 25 

Cyklotol 77/23 

1.754 

6.034 

0.0992 

4.30 

1. 10 

0.35 

H3€C 

1.891 

7.783 

0.0707 

4, 20 

1.00 

0, 30 

LX-01 

1,230 

3. 110 

0.0476 

4.50 

1.00 

0.35 

LX-04--01 

1.865 

8.364 

0.1296 

4.62 

1.25 

0, 42 

LX-07 

1.865 

8.481 

0,1710 

4.58 

1.25 

0.40 

LX-09-l 

1.840 

8.481 

0.1710 

4.56 

1.25 

0, 40 

LX-10-1 

1.865 

8.807 

0.1836 

4.62 

1.32 

0.38 

LX-11 

1.875 

7.791 

0.1067 

4.50 

1. 15 

0,30 

LX-14-0 

1.835 

8.261 

0. 1724 

4,55 

1,32 

0. 38 

LX-17-0 

1,900 

4.460 

0,1339 

3. 85 

1.03 

0, 46 

Oktol 78/22 

1.821 

7.486 

0.1336 

4. 50 

1. 20 

0.38 

PBX-9010 

1.787 

5.814 

0.0680 

4. 10 

1.00 

0.35 

PBX-9011 

1.777 

6.347 

0.0800 

4. 20 

1.00 

0, 30 

PBX-9404-3 

1.840 

8. 524 

0.1802 

4.55 

1,30 

0.38 

PBX-9407 

1.600 

5, 732 

0.1464 

4,60 

1.40 

0.32 

Pentoitt 

1.700 

5, 410 

0,0937 

4. 50 

1. 10 

0.35 

PETK 

0. 880 

3.486 

0. 1129 

7.00 

2.00 

0.24 

PETK 

1. 260 

5.731 

0.2016 

6, 00 

1.60 

0.28 

PETN 

1.500 

6.253 

0.2329 

5.25 

1.60 

0,28 

PĘTU 

1. 770 

6, 170 

0.1693 

4, 40 

1,20 

0. 25 

Tetryl 

1.730 

6.868 

0.1067 

4.40 

1.20 

0.28 

TNT 

1,630 

3.712 

0.0323 

4.15 

0.95 

0. 30 

Iłltrometan 

1. 128 

2. 092 

0.0569 

4. 40 

1,20 

0.30 
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punktowi Chapmana"Jougueta <C~J>. 2 reguły przyjmuje sie* że I 
kształt frontu fali detonacyjnej Jest płaski, sferyczny, I 
toroidalny lub stanowi ich kombinacją. Rozwiązanie konstruuje | 
się więc tylko dla samych produktów detonacji, między czołem | 
fali detonacyjnej a brzegiem. ] 
Oczywiście, takie sformułowanie zagadnienia, wymaga znajomości 1 
parametrów w punkcie C-J dla różnych materiałów wybuchowych. V '| 
tabeli nr 2 zebrano je dla tych samych materiałów, których | 
dotyczy tabela nr 1. Dane zawarte w tabeli nr 2 pozwalają | 
obliczyć gęstość 1 prędkość normalną do czoła fali 3 
detonacyjnej uTh w punkcie C“J z prostych zależności: 

<2.9X1 
< 2 . 10)1 


Równania <2,1) do <2,4) i C2.8> oraz wartości podane w tabelachJ 
112 pozwalają opisać w pełni ruch i parametry produktów I 

detonacji podstawowych, skondensowanych materiałów wybuchowych. | 

2.2. Teoria sprężysto/lepko-plastycznoścl w zastosowaniu do | 
modelowania deformacji i napędzania metalowych wkładek | 
kumulacyjnych. 1 

V pracy niniejszej zostanie omówiona aktualna i zarazem | 
dość ogólna wersja modelu opisującego zachowanie się wkładki | 
kumulacyjnej w procesie JeJ deformacji i napędzania. Pozwoliła j 
ona na osiągnięcie zgodności wyników teoretycznych z I 

eksp^^rymentalnymi w takim stopniu, iż można uznać eksperyment I 
komputerowy za praktycznie równoważny eksperymentowi \ 

fizycznemu,. Versja ta bazuje na modelu ciała sprężysto/3.epko— ^ 
pla yi_ 2 nego typu Hohenemsera-Pragera. Odpowiedni układ równań | 
różniczkowych wyrażający prawa zachowania i związki | 
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konstytutywna dla tego modelu ma następującą postać 
C93.94,95,96,97,98,99] : 

^ o diV Sr s 0 
dt 


ć ^ ^ 4. .Utl* 4- ■Sy-r-S.fy. ^2 12) 

I dt " 3r a^r ^ &Z. ^ ^ <.<i.i.^) 

o 43^ = - 4- H- + -4gg. . (2 13) 

^ +S(ftf 7^ +ć>aa^ +5^2,(^ 4^)(2.14) 

^ - id>vS?) + {|t-|^}5.z -^ęiSrr 

= 2>*- { ■^ ~ 1 divw)-,-^ ^ Slf,^ <2.16) 

^‘ = ^>^(11“ i d'-''^)- (|^~|^)S^Z->f <2.17) 

=/^f'l^'^‘lT)~ -§7)(Sr'r “ Szz.)--^^5tz. <2. 18) 

C S<j' dla \I^Y^ilSifs^ 




dz A 


ZZ (2.17) 


rz. (2.ia> 


^cl2.it ^ := 


dla sĘyy\fSil S<i 


S.i Sii-Srr^ + ^ S*»^ + 2.6^/ 


Oznaczenia: 

»3ł'| - składov/e dewiatora tensora naprężeń C Si\ ^ ** t' <S 4 .‘j ) . 

Y “ dynaml-zna granica plastycznego płynięcia, 
yU - moduł .inania, 

- współczynnik lepkości. 

Pozostałe oznaczenia sa analogiczne jak w punkcie 2,1. Równanie 
stanu uzupełniające ten model przyjęto w następuj-ącej postaci 
C1003 Ctzw. model EOS): 


p= k,K 


< 2 . 21 ) 




X - 1 f o dla x<0 ; 

a oznacza gęstość fazy ciałostałowej, która będzie 

dokładnie zdsfInloyreina w punkcie 2<3. 
yspół czynni ki ^ i. > ♦ S® i kilkunastu metali 

zebrano i przedstawiono w tabeli nr.3. 


Tabela 3. Współczynniki równań stanu dla modelu EOS. 



Metal k^C Mbar] k^C Mbar3 k^Mbar] gCg/cm^) 


Aluminium 

0.7906 

1,3250 

2.1300 

2.78 

2, 00 

■ 2łoto 

1.7970 

2.9810 

4.9320 

19.24 

2, 97 

Beryl 

1.1840 

1.9750 

2,9440 

1.85 

1, 16 

Miedź 

1.3860 

2.7490 

5,1130 

8.93 

1,99 

Magnez 

0.3511 

0.. 6376 

1.0510 

1.74 

1.42 

■ Niob 

1.6910 

2.8390 

4,3900 

8.59 

1.47 

Elkieł 

1.8790 

3.5830 

6.4300 

8.8? 

.1.93 

Ołów 

0,4774 

0,7329 

1,1220 

11.35 

2.77 

Platyna 

2.7730 

5,2350 

9.5560 

21.42 

2, 40 

Stal 

1.6480 

3,1240 

5.6490 

7.90 

2. 17 

Tantal 

1,9410 

3.1100 

4.6690 

16.65 

1,60 

Tytan 

1.2340 

1.2200 

1.1460 

4. 53 

1.09 

Uran 

1.1720 

2,3660 

4, 535 

18.95 

2.20 

Wolfram 

3.1210 

5.3180 

8,3800 

19.22 

1.54 


Model EOS pozwala również wyznaczyć temperaturę metalu z 
zależności; 



gdzie: 
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uy^zgl^dnia zależności i /i od deformacji plastycznej, 

temperatury i gęstości i 100] . Zależności te 
^^^^g^j-uowano przy prędkościach deformacji charakterystycznych 
dla procesów oddziaływania fal detonacyjnych na ciała stałe, a 
%ch postaó Jest następująca; 

bp(-|2j'^®-- h{T" 500)] <2 25> 

y^(*f'ł- ^ Ymax < 2 , 26) 

y«0 T > T nn <2.27) 

yUo j^i ^ ^ 500)j < 2 . 28 ) 

Tm ' exp [22i;(l- ||)] (2.29> 


df 




gdzie; ^ - intensywność deformacji plastycznej, 

Wartości odpowiednich współczynników dla metali 
wyszczególnionych w tabelach nr 3 i 4 podano w tabelach nr 5 l 
6 , 

V chwili obecnej nie można znaleźć w literaturze równie 
ogólnych i dobrze udokumentowanych zależności dla współczynnika 
lepkości . Dlatego też zaproponowano, aby w badaniach 

własnych uwzględnić jedynie najsilniejszą zależność 
odpowiedniego czasu relaksacji T-^ od temperatury C293 . 

Wówczas współczynnik lepkości będzie dany wyrażeniem; 




F(T) - h(T- 500)] <2.32) 


Wartość 1^0 oraz postać funkcji F(T) dobi<:.Vano na podstawie 
własnych badań teoretyczno-eksperyraentalnych wykorzystując v 
tym celu zjawisko ^‘odwrotnej kumulacji*'. 
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Okazało s;l^, ±e naj lepsze wyniki uzyskana przyjmując 
postać funkcji FCT) analogiosną do te.J, która opisuje lepkość 
płynnych metali. Do podobnych y^niosków dossli autorzy pracy 
E 101 ], badając zachowanie się ośrodków porowatych. 

Tabela 5. Współczynniki vfystąpujące w modelu Steinberga- 
Guinana^ 


Metal 

Yo-105> 



Tm© 


t Mbar] 

C Mbar] 

i Mbar] 

CK] 

Aluminium 

2 . Q 0 

0.03 

2,70 

1220 

Złoto 

0.20 

0.23 

2.80 

1970 

Beryl 

3, 30 

1.20 

15.10 

1820 

Mi edż 

1.20 

0 . 60 

4.77 

1790 

Kagnez 

1.70 

0 , 50 

1.65 

1570 

Ni Ob 

a. 00 

1.40 

3.77 

1760 

Nikiel 

i, 40 

1.20 

3.55 

1950 

Ołów 

0 . 08 

0 . 10 

0.36 

2740 

Platyna 

0.30 

0.30 

6.37 

2870 

Stal 

3. 40 

2 . 00 

7.70 

1930 

Tantal 

7. 70 

1 . 10 

6.90 

1740 

Tytan 

7.10 

1.50 

4. 34 

1230 

Uran 

8 . 00 

1, 70 

3, 44 

2420 

Volf ram 

22 . 00 

4, 00 

16, 00 

1670 

Tak wlec, 

dla miedz: 

lanych wkładek 

kumu 1 acyj nych, 


przyjmowano: 

F (T) = exp 
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1o = 5-10-^ [f] 


C2.33 


Tabela 6. Współczynniki występujące w modelu Steinberga- 
Guinana,cd. 


Metal 

P 

a 

b[Mbar"^J 

hlikl' 

Aluminium 

125 

0 . 10 

8 

6 . 2 

; Złoto 

49 

0 .39 

4 

3.2 

Beryl 

81 

0.22 

2 

2.6 

Miedź 

36 

0. 45 

3 

3,8 

Magnez 

7000 * 

0 . 10 

10 

4 .8 

Niob 

5 

0.20 

1 .4 

0 . 0 

łTlkiel 

46 

0.53 

2 

3.4 

Ołów 

110 

0,52 

14 

12,0 

Platyną 

20000 

0.0002 

3 

1.4 

Stal 

40 

0..35 

3 

4.5 

Tantal 

10 

0 . 10 

2 

1.3 

Tytan 

7S0 

0. 065 

i 

6.2 

Uran 

2. 7 

0.26 

1 

8 , 1 

Wolf ram 

7.7 

0 . 13 

1 

1.4 


Dla pozostałych metali wzór typu <2.33) musi być 
konstruowany na drodze dodatkowych, badari. teoretyczno-* 
eksperymentalnych. 

Przytoczone w tym punkcie równania oraz .opis własności 
materiałowych zamykają problem matematyczno-fizycznego opisu 
dynamicznych deformacji metali obciążanych materiałami 
wybuchowymi, w obszarach, w których pozostają one litymi 
ciałami stałymi lub cieczami, Poprawny jakościowo i ilościowo 
opis obezarow, w których tworzą się szczeliny(pustki>, wymaga 
jeszcze dusisz ej rozbudowy tego modelu, czemu poświęcony będzie 
następny punkt 2.3. 
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—. . ■i.ijdel two*^ si^ i wzrostu szczelin w deformowanycli 

i adka ch ku luu 1 ac y j nych, 


Wyniki badań eksperyowntalnych wskazują, że szereg zjawis 
towarzyszących procssOB kumulacji będzie bardzo trudno opisad 
ramach klasycznej mechaniki ośrodków ciągłych. 

Można do nich zaliczyć: rozdzielanie słą matarlału wkładt 
na strumień 1 zbitek, odrywanie się strumienia od zbitka, 
fragmentacja strumienia, fragmentacja brzegów wkładek 
kumulacyjnych, rozrywanie się nieoptymalnie zaprojektowanych 
wybuchowo kształtowanych pocisków. Itp. Adekwatny do 
rzeczywistości opis matematyczno-fizyczny musi stwarzać 
możliwości mqdelov/ania tych efektów, które wiaża się z 
powstawaniem mikro jak i makro nieciągłości ośrodka. 
Hleclagłoścl te przyjęto nazywać szczelinami gdy tworzą one 
wydłużone magistrale pęknięć lub pustkami <poraml>, gdy ich 
geometryczne kształty są bardziej regularne, np. zbliżone do 
Sfer lub elipsoid. 

Zdecydowana wl^kszo^ć prac dotycząca tej probleiaatyki 
koncentruje sle wokół sagadnleA powstawania i 

rozprzestrzeniania sIę poJedyiiczych szczelin. Z reguły traktuj^ 
SU Je Jak wolne od iiapreże:*i powierzchnie w odpowiednim 
problemie brzegowym fizyki matematycznej. Warunki progowe na | 
powstawanie 1 wzrost szczelin określa się zwykle przez podanie f 
^ ęnergil, intensywności naprężeń lub 

odk=2tał..eń. Podejście takie ma już w chwili obecnej dość 
bogatą literaturę <np. C 102.i-108) 1 ,aie niestety nie Jest ona ^ 
zbyt pomocną przy komputerowym modelowaniu złożonych, 

niewyldealLzowanych problemów. Wynika to z dwóch zasadniczych ^ 
powodów. 

Po pierwsze, śledzenie rozwoju pojedyftczej szczeliny i 
o.ipow..ednie JeJ modelowanie wymaga skomplikowanych zabiegów 
rozdzielania siatki numerycznej i tworzenia nowych brzegów 
wewnątrz rozważanego obszaru. Przedsięwzięcia takie stają się ' 
szybko nierealne technicznie Jeśli liczba niezależnych szczelin 


y^zrastaó Lub też kierunki ich rozprzestrsreniania si^ 
niezgodne z układem wązłów sieci numerycznej. 

Po drugie* w świetle nowszych badah C109,1101 zagadnienia 
piszczenia struktury plastycznych metali (miedź, żelazo armco, 
nie można sprowadzić tylko da analizy propagacji 
jaakro^sczelln. Proces zniszczenia zapoczątkowują bowiem licznie 
powstające, niezależne mikropąknlącia różnych kształów i 
rozmiarów,chaotycznie zapełniające niszczony obszar. Dopiero w 
pdżnfedszych fazach zjawiska, na skutek łączenia sie 
mikropekniód, powstają makroszczellny. Bardzo istotny Jest przy 
fakt, że powstające mikropekni^cia mogą silnie zmieniać 

wytrzymałościowe ośrodka, w dość znacznych obszarach. 

Dlatego też przy modelowaniu komputerowym 
wielowymiarowych, złożonych procesów z uwzględnieniem 
niszczenia struktury metali, podejście polegające na śledzeniu 
rozwoju pojedyAczych makroszczelin Jest rzadko stosowane 1321 i 
uznać Je za nleperspektywtczne z wymienionych wyżej 
technicznych t fizycznych powodów. 

V chwili obecnej zaczyna dominować w taj dziedzinie inny 
sposób opisu powstawania i wzrostu objętości zarówno mikro Jak 
I inakroszczelin. Kówląc ogólnie, polega on na tym, że baz 
wnikania w strukturę pojedyAczych szczelin ,ocenia sle jedynie 
Ich średnią objętość w danym elemencie ośrodka i w zależności 
od niej modyfikuje mechaniczne właściwości metali. 

Najprostszy model tego typu przedstawiono w pracy Clili. 
Sprowadza się on do formalnego ograniczenia; ciśnienia, granicy 
plastycznego płynięcia i modułu' ścinania w obszarach, w których 
gęstość spada poniżej pewnej krytycznej wartości. Podejście 
"tafcie ma charakter statyczny 1 zapewnia jedynie, że zmiany tych 
wielkości odbywają się w fizycznie uzasadnionych granicach, gdy 
w ośrodku zaczynają powstawać szczeliny. Ze względu na swą 
prostotę sposób ten był wykorzystywany w początkowych etapach 
niniejszej pracy C31 . 

Bardziej zaawansowane, dynamiczne modele przedstawiono w 
pracach (42-5-46,92, 101, 109, 110,1123 . Niestety, prace te 
pojawiły się stosunkowo niedawno i nie doszła jeszcze do 
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określonego ujednolicenia poglądów w tej probleinatyce. Równania 
opisujące dynamikę wzrostu objętości szczelin oraz icb wpływ na 
własności wytrzymałościowe metali różnią się u poszczególnych 
autoi ów zarówno formą Jak 1 wartościami liczbowymi odpowiednich 
współczynników, 

V tej sytuacji zdecydowano się na zaadoptowanie podejścia 
przedstawionego w pracach E44,45}, ze względu na jego dość 
szeroką 1 zakończoną pozytywnymi rezultatami weryfikacją 
eksperymentalną <dla żelaza armco i aluminium). Okazało się 
później, że wybór ten był dość trafny, a odpowiednie 
współczynniki liczbowe dla miedzi znaleziono w oparciu o własne 
prace teoretyczno-eksperymentalne. 

Układ równań opisujący dynamikę powstawania i wzrostu 
objętości szczelin przyjęto więc w następującej postaci; 


^ = - k (p)[|p| - ei, 1 (Vc ^ Vco) (2.34> 

•^=0 dl, Ip|< ,2.35) 


Znajomość objętości szczelin pozwala łatwo wyliczyć-gęstość 
fazy clałostałowej ośrodka z relacji: 

1 = V + J_ 

S §S (2.3f 

gdzie; 'Vi, - objątośi szczelin w Jednostce masy ośrodka, 

5 - gęstość średnia, 

“ gęstość fazy ciałostałowej. 

Wartości stałych współczynników dla miedzi wynoszą: 




Vc;«5-<0‘5 


stała ma sens fizyczny początkowej objętości szczelin w 

niezaburzonym metalu, ■ 
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Yyraźenle określa wartość progową średniego naprężenia, 

pQ przekroczeniu którego objętość szczelin rośnie < 0 > lub 

maleje < p >0 ). 

Ze związku <2,36) wynika, że w obszarach , w których 
powstały szczeliny < 0 >, ośrodek ma strukturę dwufazową - 

jego objętość właściwa jest sumą objętości właściwej szczelin 1 
ciała stałego, V takim ośrodku równanie stanu 1 model 
Steinberga mogą opisywać tylko stan fazy clałostałowej . Dlatego 
też ciśnienie oraz wielkości takie, jak; Y ^ t należy 
^uzależnić tylko od gęstośćt fazy clałostałowej §5 , a nie od 
gęstości średniej g (uwzględniono to już w zależnościach 
podanych w punkcie 2 , 2 ), 

Jak już wcześniej wspomniana pojawiające się szczeliny 
modyfikują własności wytrzymałościowe ośrodka. W modelu 
matematyczno-fizycznym można ten efekt uwzględnić wprowadzając 
odpowiednie ograniczenia wielkości Y , yu. i . zależność! od 
objętości szczelin. Ośrodek ze szczelinami będą wówczas 
ćharakteryzować efektywne wielkości Y^t yU."*" 1 wyliczane 

następująco; 

y^»r-Gac}; ( 2 . 37 ) 

Dla miedzi przyjęto; 

G{Vc) = exp(-?S'Vd) . ^2,38) 

Postać funkcji G(v; określono w oparciu o własne badania 
teoretyczno-eksperymentalne, gdyż żadna z literaturowych 
postaci &(\^)nie dawała dostatecznie zadowalających rezultatów. 
Stosunkowo najbliższe propozycji określonej wzorem <2.36), są 
funkcje &fVc) podane przez autorów prac C42] i t45]. 

Znajomość objętości szczelin w poszczególnych elementach 
ośrodka jest również bardzo pomocną przy formułowaniu kryteriów 
fragmentacji ciał na niezależne części* Do problemu tego 
wrócimy w rozdziale 4, gdyż łatwiej Jest go dyskutować na bazie 
rozwiązań konkretnych problemów. 
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2.4. Warunki początkowo-brzegowe. 


Warunki początkowe opisu łące rilezabu rżany stan wkładki I 
kuitiulacylnej 3prav/adi:aJą dc przyjęcia w chwili początkowej 

następujących wartości zmiennych zależnych: ^ 

5(/r^2,0)^go ; >u(or^2!L, 0) = 0 ; 

; S^iKl^nZjOl^O J Sitc(r^z,o)^Oj 

o i 

Przy pomocy kodu komputerowego, który bidzie omówiony w | 
następnym ro<.d 2 iale, można rozwiązywać szeroką gamę różnych 
problemów brzegowych. Dlatego też warunki brzegowe należy | 
sformułować dość ogólnie, abstrahując od konkretnych | 

przykładów. | 

Tak wiec n.ą wszystkich swobodnych powierzchniach wkładki I 
, przyj mujemy: 

6;('r,z.t)l =0 ; 6V('r',z.łr)| =0 „ 1 

lł<K*.+ )=0 > ^ <2.4a| 

I 

gd^le: 6'^ 1 Sx oznaczają odpowiednio naprężenia normalne 1 

styczne do powierzchni wkładki, | 

I 

Natomiast na powierzchniach swobodnych produktów detonacj| 
“ "fż. mamy; | 

Kd powierzchni kontaktu; produkty detonacji - wkładka 
kumulacyjna, zakładano warunek tzw. swobodnego poślizgu; 

-o ] 

\ . <3..» 

gdzie; (n‘,ait)'0- równanie powierzchni kontaktowej, 

■ hormalne do powierzchni kontaktu, składov/e J 
prędkości masowych wkładki i produktów 
detonacji. I 



Zgadnie z przyjętą metodą opisu procesów detonacji MV, 
ę. propaguj ący sie front fali detonacyjnej Jest modelowany 
powierzchnią silnej nieciągłości o zadanym kształcie 1 
prędkości przemieszczania się. Na powierzchni tej zakłada się, 
że odpowiednie zmienne zależne przyjmują wartości odpowiadające 
punktowi C”J; 

p~p+i > C2.43> 

Wielkości i ^ 4 ^ należy w każdym punkcie dobierać tak, aby 

wypadkowy wektor prędkości (U^^ltę^^był skierowany normalnie do 
powierzchni frontu. 

Warunki początkowo—brzegowe zamykają sformułowany w 
punktach 2.1 2,3 matematyczno-fizyczny opis problemu. 
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Rozdział 


STRUKTURA I PODSTAWOyS HLEMEHTY KODU NUMERYCZNEGO 
DO ANALIZY SPRZĘŻONYCH PROCESÓW DETONACJI»DBFORMACJI 
I NAPĘDZANIA CIAŁ. 


Złożoność sforinułowanych w poprzednim rozdziale równań | 
problemu sprawia> ź© nie dają one rozwiązywać metodami | 

analitycznymi. Szansą na uzysKanl© przybliżonych rozwiązań ty| 
równań stwarzają Jedynie metody fizyki komputerowej, do któryl 
zalicza sl© również metoda, która bądzle omówiona w ninlejszyl 
rozdziale. Szersze omówienie tej metody jest podyktowane tym, ; 
że nie ma ona swoich bezpośrednich odpowiedników J 

literaturowych, a rozwiązania poszczególnych problemów mają 
własny, oryginalny charakter. 

Przy analizach numerycznych złożonych procesów fizycznych 
(zwanych często symulacjami komputerowymi) osobnym problemem | 
staje sią ocena poprawności otrzymywanych wyników. ¥ naszym | 
przypadku, v^ielowymiarowość zagadnienia oraz nieliniowość 
układów równań powodują, że udowodnienie analitycznymi metadaij 
stabilności lub zbieżności schematu różnicowego Jest 
praktycznie niemożliwe. Jedynym, praktycznym sposobem oceny 
poprąvmoścl wyników, Jest w takim przypadku konfrontacja ' 

wyników symulacji komputerowej ze znanymi rozwiązaniami 
niektórych szczegółowych problemów lub wprost z wynikami 
eksperymentów. Konfrontacje tego typu pokazały, że proponowan 
metoda oddaje poprawnie zarówno Jakościowe jak i ilościowe f 
charakterystyki wszystkich dotychczas analizowanych zjawisk. I 
Przejdźmy wląc do JeJ bezpośredniego scharakteryzowania od | 
strony ideowej i częściowo realizacyjnej. I 

3.1. Idea proponowanej metody numerycznej, I 


Wspólną cechą typowych algorytmów wykorzystujących 
współrzędne Lagrange’a lub Eulera Jest to, że tworzona tam si^ 
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dzieli ośrodek na poszczególne komórki, w których 
-S^^^enrowadza si© odpowlfcanie bilanse masy, p©du i energii. Na 


podstawie określa s±^ następni© wartości poszukiwanych 

metoda różni sl© od tych metod dość zasadniczo 
^ etapie tworzenia samego pojąc la sieci numerycznej . 
^^|^ 2 ada si© w niej bowiem, że specyficzną st©ć numeryczną 
wybrane w chwili początkowej punkty 
^^gaterialne, przemieszczające si© wraz z ośrodkiem. Rozwiązanie 
ttiUioeryczne sprowadza slą w takim przypadku do obliczania w 
4^^^.^jnych chwilach czasu położeń 1 wartości zmiennych zależnych 
punktach. Tak wląc, bilanse odpowiednich wielkości w 
poszczególnych “komórkach", zastępuj© sl© w tej metodzie 
J^hliczantem trajektorii i parametrów punktów materialnych. 
^Vunkty te w proponowanej metodzie nie. tworzą trwałego 
-sąsiedztwa, tak Jak to ma miejsce na przykład dla wązłów sieci 
tagrang 0 *a. Mogąc dowolnie zmieniać położenia wzglądem siebie, 
punkty te tworzą specyficzną, zmieniającą sl© w czasie 1 
przestrzeni sieć z "wymiennymi" węzłami. Ta ”wymlenność'* wązłów 
-pozwala na modelowanie dowolnie dużych deformacji ośrodka.’ Jest 
tp podstawowa zaleta tej metody, którą można by opisowo 
określić jako; metod© Lagrange'a przystosowaną do obliczeń w 
warunkach bardzo dużych deformacji. 

Przejdźmy zatem do omówienia konkretnej realizacji takiej 
"i^tody. Załóżmy, że badane ciało zostało pokryte siecią punktów 
materialnych, w których znane są wszystkie wartości zmiennych 
zależnych w chwili czasu ¥ybierzmy Jeden z takich punktów 1 

:: przyporządkuj my mu par© liczb całkowitych <L,K>, która będzie 
go Identyfikować. Liczby L i K można utożsamić z 
:>:lągr a nge'o wsk i mi współrzędnymi danego punktu. Zgadnie z tymi 
^GZnaczaniaml położenie punktu <L,K) na płaszczyźnie <r,z> w 
chwili t*' opisuje para liczb <r ,z ), a wartości 
zmiennych zależnych w tym punkcie opisują odpowiednio wielkości 
typu -<F i ) . Położenia punktów tworzących lokalne 

sąsiedztwo punktu <L,K> oznaczymy następującej Cr JJ, tZ m 
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m Wskaźnik ‘♦m'* numeruj© punkty sąsiednie* a ”N'* określa 

ich liczbę. Wartości zraiannych zależnych w punktach sąsiednich 
można analogicznie oznaczyć -<F£ . Elementarny i zarazem 

podstawowy problem sprowadza sie teraz do odpowiedzi na 
pytanie: w jaki sposób po upływie czasu <t -t*^= A t> 

obliczyć połażenie i nowe wartości zmiennych zależnych w 
punkcie (L*K> ?, 2auważmy w tym celu* że równania problemu 
można przedstawić w następującej ogólnej postaci: 

Gi (kk.Fj ,-f^^) (3.1 

gdzie: M; Z 

^ T ) X2."Z J 

a M Jest liczbą zmiennych zależnych. Pochodna substancjalna 
opisuje Jak wiadomo zmianą parametrów w poruszającej sią 
cząstce materialnej, a wiąc forma zapisu równań typu <3,1> 
sugeruje Już cząśclową odpowiedź na powyższe pytanie. Można 


bowiem wartość funkcji F, 
t wyliczyć następująco: 


punkcie i <L,K> w chwili 


:IFVJI« . (GO", 


Postać prawej strony róvmania <3,2> świadczy o tym, że 
zdecydowano się na wybór prostego,Jawnego schematu 
numerycznego. Ma początku pracy wykorzystywano również bardziej 
złożone, wielofcrokowe typy schematów C33. Okazało slą Jednak, 
że im bardziej złożone zagadnienia trzeba rozwiązywać, tym 
coraz bardziej problematyczne staje się korzystanie z 
dokładnieJszych** schematów. Z jednej bowiem strony, w 
złażonych, wielowymiarowych problemach symulacyjnych zatraca 
sle możliwość subtelnej oceny dokładności schematów, a z 
drugiej strony owe *'dokładniejsze** schematy mogą zbytnio 
komplikować cały algorytm i czynić go mniej efektywnym. 
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Aby skorzystać z równania <3.2> musimy umieć wyliczyć 
llinkcją ( i ^ Jest to problem taki prosty, gdyż w 

punkcie <L,K) nie znamy gradientów , Dysponujemy Jedynie 

ihformacjami o wartościach funkcji F^‘ w dyskretnych, 

?iieregulamie rozłożonych wokół punktu <L,K> punktach 
Sąsiednich < IT, .Z. ^ ^ . ¥ pracy niniejszej zaproponowano 
następujący sposób obliczania gradientów który można 

prześledzić posługując sle rysunkiem nr 8. * 


Rys.0. Przykładowa ilustracja położeń: punktu obliczanego 
L|Ic> t punktów sąsiednich < T JJ, t 21 Si > i 
punktów sąsiednich interpolowanych <T^rn,21 fti >. 

Parametry 1 położenia nieregularnie rozłożonych punktów 
sąsiednich interpolujemy liniowo na okrąg o promieniu R, 
którego środ,kiem Jest punkt <L,K>. Położenia i parametry w 
punktach Interpolowanych można łatwo wyznaczyć: 


(tw - T (2iJ« 
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<3. 6) 


Promień R wygodnie Jest dobierać tak aby był on bliski 
średniej odległości miądzy punktem CL,K)t a punktami 
sąsiednimi. Jest on w tym schemacie niejako równoważny 
klasycznemu krokowi przestrzennemu. Interpolacja C3.3) -5- (3.6) 
ma na celu nadanie wszystkim punktom sąsiednim tej samej ”wagi" 
względem punktu obliczanego. 

Przyjmijmy dalej, że funkcją można wokół punktu CL,K> 

lokalnie zaaproksymować funkcją liniową: 

= (Fj)”,k. +a(r-'ru,K) ^ bCz-“2.",iŁ) (3.7> 

Współczynniki alb można dobrać tak, aby płaszczyzna <3.7> 
możliwie najlepiej aproksyioowała znany, dyskretny rozkład 
wartości funkcji ( * 

Utwórzmy w tym celu wyrażenie: 

A(a.b)=i^ ..3.8> 

Wyrażenie to Jest sumą kwadratów różnic mlądzy wartościami 
fi'nkcjl C w punktach ^ wartościami Jakie w 

tych punktach wynikają z wyrażenia (3.7>, Aby różnice te były 
najmniejsze, musi zachodzić: 

<3.9). 

<3.10> 

Wprowadzając następujące oznaczenia: 
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Z (3.7> wynika bezpośrednia, że poszukiwane pochodne 
s t rz e n ne są r ó wne c 


(3.11> 


<3.12) 


<3,13) 


<3,14) 


<3.15) 


<3.16> 


<3,16) 


Przedstawiony sposób oceny gradientów wokół punktu 
. obliczanego jest niezawodny w, praktyce numerycznej, tzn. nie 
i.daje nigdy nie fizycznych wielkości gradientów. Jedynym 
waru jaki musi być przy tym spełniony jest odpowiednia 

gęstość kątowych rozkładów punktów sąsiednich, co bądzie 
szerzej omówione w następnych punktach. 

Z <3.2) możemy Już teraz bezpośrednio wyliczyć wszystkie 
interesujące nas wartości zmiennych zależnych w| punkcie <L,K) w 
"chwili 

B A+[GiK,Fi,.|Jt)K.. -1.. 


ITowe położenie punktu <L,K> wyznacza się następująco: 
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r - r E,k u • M 
Z.1*K ~ 2 , 1 , K <- O-C.K ■ At 


Okazało sl^ jednakże, że tak pomyślany schemat numeryczny ^ 
Jest niestabilny, tzn. generuje on szybko narastające j 

nleflzyczne fluktuacje w wyniku czego szybko rozpada sie całe ^ 
rozwiązanie. Fluktuacje te Jest bardzo trudno wyeliminować tak ^ 
aby równocześnie nie dochodziło do zbytniego, numerycznego ź 

wygładzania ^rozwiązań. W wyniku wielu badah znaleziono = 

ostatecznie postać tzw. ”zewnętrznej ,nielokalnej dyfuzji 
numerycznej” , która stabilizuje rozwiązanie i nie rozmywa 
frontów falowych w stopniu większym niż, np.pseudolepkość♦ 
Zagadnieniu temu poświęcony będzie również osobny punkt 
niniejszego rozdziału. W tym miejscu można jedynie nadmienić, 
że formalnie problem ten sprowadza sle do odpowiedniej 
modyfikacji schematu (3*19) zastosowanego do równah zachowania:^ 

pędu. Ą 

Jeżeli przedstawioną powyżej metodę zastosować kolejno do 
wszystkich punktów materialnych, pokrywających badane ciało 1 
powtarzać Ją w kolejnych chwilach czasu, to dostaniemy 
poszukiwane rozv^iazanie problemu w dyskretnej , czasowo- ^ 

przestrzennej formie. ' .'M 

3.2. Wybór punktów sąsiednich. I, 

Z przeprowadzonych dotychczas rozważań wynika, że w | 

trakcie obliczeń w każdym punkcie musimy dysponować informacją :| 
o aktualnie najbliższych względem niego punktach sąsiednich 1 . :| 
Ich parametrach. Odpowiednia obróbka tej informacji pozwala | 
określić niezbędne gradienty pól prędkości i naprężeń. I 

V tej sytuacji należy w pierv/szej kolejności udzielić | 

odpowiedzi na następujące pytanie: z ilu 1 w Jaki sposób I 

wybieranych punktów należy tworzyć zbiór punktów sąsiednich dlai 
danego punktu ?. Na pytanie to nie można niestety udzielić i 
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, jednoznacznej odpowiedzi C491, Hożna jedynie odwołując 
^lę do określonych analogii z klasyczną siecią numeryczną, 
^/skazać ,na dość dobre i sprawdzone w praktyce rozwiązanie tego 
robleiHU. Pochodne przestrzenne w poszczególnych węzłach 
j^ggycznej sieci numerycznej wylicza się w oparciu o informację 
2 ^Wartą w czterech węzłach sąsiednich. Kątowy rozkład tych 
względem węzła, w którym aktualnie prowadzone są 
Obliczenia charakteryzuje kąt = 90®. Indeksy 

J ''> oznaczają numery dwóch, kolejnych w kątowym rozkładzie 
vfęzłów^, tak jak pokazano to przykładowo na rys. 9. V omawianej 
metodzie rozkłady kątowe punktów sąsiednich będą nieregularne, 
ale możemy zarżądać, aby ilość 1 sposób ich wyboru zapewniał 
^^pełnienle warunku Cfi;/ ^ 90®. 



Sys,9. Ilustracja metody wyboru punktów sąsiednich 


l Z rys.9 widać, że warunek ten można spełnić jeśli z punktu 

^ <r?,K ♦2ŁU|ic> podzielić płaszczyznę <r,z> na osiem równych 

[ sektorów kątowych l w każdym sektorze wybrać punkt najbliższy 

j ^względem punktu <L,K), Wówczas, nawet przynajmniej korzystnym 
I rozkładzie kątowym <p. 112 na rys. 9) będzie zachodzić (^^<90®, 
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Vybór oiśmlu punktów do tworajenia zbioru najbli^^jszych 
sąsiadów dla danego punktu w oparciu o wyżsj sforzrsułoyrane 
kątowo-odległości owe kryterium można wi^jc w pevmyin sensie uzna^ 
za "optymalny”. Mniejsza liczba punktów sąsiednich zwiększa 
bowiem niebezpieczeństwo pojawienia sie tzw. niestabilnośól 
kątowych ( ^90®) . Natomiast zwiększenie liczby punktów 

sąsiednich powyżej 6 prowadzi jedynie do niepożądanego 
zwiększania obszaru wzajemnych oddziaływań między punktami. 

Ze wzglądu na dużą czasochłonność operacji ustalania 
.najbliższego. sąsiedztwa, nie można jej przeprowadzać na każdymi 
kroku czasowym. V zależności od szybkości 1 charakteru 
deformacji ośrodka, odstąp czasowy między kolejnymi operacjami! 
uaktualniania zbiorów najbliższych sąsiadów, powinien być 
dobierany tak, aby przeciętnie wymianie podlegał nie więcej nlżi 
jeden spośród ośmiu punktów,'Jest to podyktowane dążeniem do 
zminimalizowania fluktuacji numerycznych, które zawsze 
wprowadza "operacja wymiany sąsiadów**. Operacja ta zapewnia 
wspomnianą Już w poprzednim rozdziale wymienność" węzłów sieci 
numerycznej, umożliwiającą prowadzenie obliczeń w warunkach 
skrajnie dużych deformacji. 

Czasochłonność obliczeń można znacznie zmniejszyć również ; 
na innej drodze, mianowicie każdemu punktowi można 
przyporządkować pewien podzbiór M punktów i tylko z'niego 
w>^blerać najbMższe osiem punktów sąsiednich, a nte przeglądać 
zbioru wszystkich punktów. Podzbiór ten można łatwo utworzyć na 
b.azie aktualnych punktów sąsiednich wziętych wraz z ich 
własnymi punktami sąsiednimi. Podzbiór M może więc w takim 
przypadku zawierać maksymalnie 72. punkty. Operacja jego 
tworzenia Jest technicznie dość prosta, jeśli w każdym punkcie 
zakodowana jest informacja o numerach Identyfikacyjnych jego 
punktów sąsiednich. Konieczność kodowania takiej Inforir^cji 
wynika również stąd, że zbiór punktów sąsiednich musi być w 
jakiś sposób zapamiętany między kolejnymi operacjami Jego 
uaktualniania. 


Dla punktów tworzących brzeg rozważanego obszaru Ca 
niekiedy 1 dla niektórych punktów wewnętrznych) powstawać 
indzie problem braku punktów sąsiednich w niektórych sektorach 
Icątowych. V takim przypadku najwygodniej podzielić proces 

zbioru punktów sąsiednich na dwa etapy. W pierwszym 
. itapie ż każdego sektora, w którym znajdują się punkty z 
podzbioru M, wybieramy punkt najbliższy względem danego punktu. 
Brakując© punkty Cdo ośmiu) można w drugim etapie dobrać z 
pozostałych punktów podzbioru H rezygnując Już z kryterium 
kątowego wyłącznie na rzecz kryterium najmniejszych odległości. 

Ba rys.9 takim dodatkowo dobranym punktem jest punkt nr 7, gdyż 
sektor nr VII jest sektorem pustym. 

Oczywiście, możliwość pojawienia sle sektorów pustych dla 
■^punktów wewn^ętrznych może prowadzić do naruszenia kryterium 

(90*=* i trzeba temu w określony sposób przeciwdziałać aby 
uniknąć pojawienia się tzw. niestabilności kątowych. Problem 
ten zostanie szerzej omówiony w następnym punkcie, a na 
'zakończenie tych rozważań należy zwróció jeszcze uwagę na pewne 
praktyczne aspekty poprawności realizacji procesu wyboru 
punktów sąsiednich. V praktyce zdarzać się bowiem będą 
pr^ że punkty spełniające formalne, kątowo-odległościowe 

lWunkl. nie będą mogły fizycznie tworzyć sąsiedztwa danego 
•punktu. Może to wynikać na przykład z rozdzielenia tych punktów 
brzegami rozważanych obszarów. Dwa tego typu przykłady pokazana 
,na rys,10. 

punkty takie, Jak np. oznaczone numerami 1 1 2 na rys. lOa) 1 lOb) 
iuszą być eliminowane ze zbioru punktów sąslednlcti wzglądem 
punktu <L,K) aby zapobiec pojawianiu slą nlellzycznycb efektów. 
Technicznie, operacje takie muszą wykonywać odpowiednie 
algorytmy śledzące usytuowanie punktów sąsiednich i linii 
brzegowych. 
















/ 
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® deionac^ 
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1 

(l.K^ próżnia 


•«2 

^ linia brzegowa 

oirodek 

/rkumulacyjńa : 



xys, 10. Punkty <1,2> ellioinowane zg zbioru punktów 
sąsiednich punktu <L»K>. 

3.3. KiestablInoścl zbliżeniowe i kątowe. 


Jedną z charakterystycznych i zarazem dość nietypowych 
cech omawianej metody numerycznej Jest to, że sieć przestrzennąl 
tworzą w niej nieregularnie rozłożone■punkty przemieszczające 
się vfraz z deformowanym ośrodkiem. Stwarza to określone 
trudności przy próbach zdefiniowania tzw. kroku przestrzennego," 
odpowiedniego dla tej metody, 2 kolei znajomość tego kroku JestJ 
bardzo ważna przy doborze kroku czasowego zapewniającego 
stabilność metody numerycznej. 

V trakcie badań testowych okazało sią , że z punktu 
widzenia warunków stabilności metody, rolą kroku przestrzennego? 
odgrywa odległość miedzy danym punktem, a najbliższym punktem 
sąsiednim. Oznacza to, że jeśli przy danym kroku czasowym, w 
pewnym miejscu rozważanego obszaru, dwa punkty zbliżą sle na 
odległość mniejszą od pewnej wartości krytycznej (do której 
dostosowano ten krok czasowy), to nastąpi w tym miejscu lokalne. 
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przekroczenie warunków stabilności. Objawia sie to szybkim 
narastaniem niefIzycznych wielkości wszystkich parametrów w 
otoczeniu tych punktów, a w dalszej kolejności prowadzi do 
rozpadu całego rozwiązania. HiestabiLneść tego typu można wiec 
na^^/ać ‘♦niestabilnością zbliżeniową”, gdyż pojawia sle ona przy 
2 v,ytnim (dla danego kroku czasowego) zbliżeniu do danego 
punktu, któregoś z Jego punktów sąsiednich. 

Oczywiście, można w takiej sytuacji próbować unikać tego 
typu niestabilności zmniejszając krok czasowy odpowiednio do 
zmniejszających sle odległości między punktami. Rozwiązanie 
takie Jest Jednak bardzo nieefektywne. Pojedyncze punkty mogą 
bowiem spowodować tak znaczne zmniejszenie kroku czasowego, że 
że względów technicznych sens straci dalsze kontynuowanie 
obitczeh. 

0 wiele bardziej racjonalny jest inny sposób 
zabezpieczenia algorytmu przed niestabllnościaml zbliżeniowymi, 
“Można bowiem ustalić pewien minimalny promień taki, że 

zbliżenie dwóch punktów na odległość mniejszą od Rmir>t powoduje 
wykluczenie jednego z tych punktów z dalszych obliczeń. 
Równocześnie z tym wykluczeniem musi oczywiście następować 
ustalenie nowych zbiorów najbliższych punktów sąsiednich. Można 
również nie wykluczając żadnego punktu z obliczeń, wykluczyć 
Jedynie oba punkty wzajemnie z© zbiorów swych najbliższych 
punktów sąsiednich. V ten sposób punkty te, których odległości 
są mniejsze od nie "widzą się” wzajemnie 1 nie 

“oddztaływują" bezpośrednio ze sobą. Oba sposoby eliminowania 
niestabilności zbllżeniow^^ch są w równym stopniu skuteczne i 
oba mogą być wykorzystywane w praktyce numerycznej. 

V trakcie omawiania sposobu wyboru punktów sąsiednich 
uwypuklono fakt, że przy nieregularnym kątowym rozkładzie 
punktów sąsiednich trzeba Jednak dbać o to, aby rozkład ten 
spełniał warunek Rozkład kątowy nie spełniający tego 

warunku będzie prowadził do określonego zafałszowania 
obliczonych na Jego podstawie pochodnych przestrzennych. Jeśli 
bowiem w dużym sektorze kątowym < >90®) nie ma punktów 
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sąsiednich, to nie nsamy taia róvmleiź informacj i o varto.^ciach | 
odpowiednich pól. Tyjnczasero metoda numeryczna formalnie ^ 

•*wyinterpoluJe« te inforinacje na podstawie danych zawartych w j 
pozostałych sektorach. Procedura taka prowadzi Jednak zawsze d^ 
stopniowej utraty - początkowo dokładności, a w późniejszych | 
fazach stabilności lokalnej i globalnej (analogiczna sytuacja 
w>'stąpiłaby dla klasycznej sieci numerycznej, gdyby centralne | 
pochodne przestrzenne, w sposób losowy zastopować lokalnie | 
pochodnymi lewo lub prawostronnymi). Ponieważ niestabilności 
związane są w omawianej metodzie z nieodpowiednim rozkładem I 
kątowym punktów sąsiednich, nazwana Je '*ntestabllnościami 1 
kątowymi". J 

Praktyka numeryczna pokazała, że kryterium $90*=* moźel 

być w określonych warunkach nieznacznie naruszone. Stopieh w 4 
Jakim można Je naruszyć, zależy jednak od położeń 1 stanów | 
poszczególnych punktów, a w szczególności od poziomu | 

"zewnętrznej dyfuzji numerycznej" . Problemu tego nie można | 
ściśle i Jednoznacznie rozstrzygnąć i dlatego w praktyce 
bezpieczniej Jest starać się świadomie nie naruszać warunku f 

Zaproponowany sposób wyboru punktów sąsiednich zasadnicza " 
gwarantuje spełnienie tego warunku. Powstają Jednak i takie - 
sytuacje, w których Jego spełnienie Jest fizycznie, niemożliwe, ^ 
a utrzymanie stabilności obliczeń wymaga specjalnego ; 

przeciwdziałania. Ilustruje te sytuacje rys.11. ^ 

Ka rys.11 a)numerami 1,2 i 3 oznaczono punkty tworzące 
brzeg rozważanego obszaru , a literami A i B przykładowe punkty 
wewnętrzne. Z punktów A i B najbliższe punkty brzegowe widać 
odpowiednio pod kątami i Dla punktu A zachodzi 

^ sytuacja taka nie wymaga żadnych interwencji. 
Hatomiast dla punktu B mamy a w sektorze tym nie ma 

Już fizycznie żadnych innych punktów. Sytuacja taka prowadzi 
nieuchronnie do pojawienia się niestabilności kątowej. Jedynym 
sposobem JeJ uniknięcia Jest usunięcie punktu B z dalszych 
obliczeń, z równoczesnym utworzeniem nowych zbiorów punktów 
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Rys,11, Przestrzenne konfiguracje punktów sąsiednich 

sytuacje sprzyjające powstawaniu niestabilności 
kątowych. 

eąsisdnlch, Aby zapobiec zbytniemu rozrzedzaniu sieci 
przestrzennej skutkiem takich zabiegów, można w miejsce 
wyeliminowanego punktu wewnętrznego B, wprowadzić między punkty 
213 nowy punkt brzegowy B*. Parametry tego punktu mogą 
wynikać z liniowej interpolacji miedzy parametrami punktów 2 1 
3. Działanie takie prowadzi do lokalnego zagęszczenia sieci 
punktów brzegowych tam, gdzie odległości między nimi stały się 
żbyt duże. Zioniejsza się w ten sposób prawdopodobieństwo 
pojawienia się w tym obszarze dalszych niestabilności kątowych 
1 jednocześnie nie zmniejsza się ogólnej liczby węzłów 
numerycznych. 

Ifatoratast na rys. II b) pokazano typowe sytuacje w Jakich 
Ifiogą znaleźć się punkty wewnętrzne A,B i C przy zbliżaniu się 
do osi symetrii. Jeżeli w sektorach kątowych skierowanych na oś 
Symetrii znajdują się dwa punkty sąsiednie, tak Jak to ma 












lalejsce w przypadku punktu A, to nie zscłiodzl ocsywlścle 
nlsbezpieczańst-rfc pojawienia sią niestabilności te,T;awych. 

Jednak punkt wewnętrzny znajdzie sie dostatecznie blisko osi, | 
to w sektorach tych ooie znaleźć się Jeden punkt sąsiedni | 

<punkt B) lub może nie być ich w ogóle (punkt C). W takiej j 

sytuacji należy w obszar sąsiedztwa punktów B i C włączyć | 

dodatkowe, specjalnie w tym celu tworzone punkty B’ i C’, zwanij 
dalej punktami fikcyjnymi. Parametry punktów fikcyjnych nie sąl 
obliczane, ale zadawane w taki sposób, aby zapewnić spełnlenle| 
odpowiednich warunków na osi symetrii <u=0, 6V*=0>. Tak wiąc | 
punkty fikcyjne mogą z jednej strony zapobiec powstawaniu | 

przyosiowycb niestabilności kątowych, a z drugiej strony mogą | 
umożliwić spełnienie odpowiednich warunków brzegowych. Problem? 
ten łączy slą wlec a problemem numerycznej aproksymacji ! 

warunków brzegowych i omówiony zostanie dokładniej w następnymi 
punkcie. 

3.4. Aproksymacja numeryczna warunków początkowo-brzegowych. 


V kodzie numerycznym wszystkie zmienne określające stan i 
każdego z punktów materialnych są zapamiętywane 1 przechowywana 
w odpowiednich , wielowymiarowych tablicach. Warunki początkowe 
określają wartości wszystkich zmiennych zależnych w chwili 
początkowej, a wlec wystarczy wartości te wprowadzić przed 
rozpoczęciem obliczeń w określone miejsca tych tablic. Problem 
ten nieco komplikuje sle przy modelowaniu detonacji materiału 
wybuchowego. V tym przypadku należy bowiem zapamiętać dodatkowo 
czas dojścia frontu fali detonacyjnej do danego punktu oraz 
umieć ro2r,-łżnlać punkty przed i za frontem tej fali. Problem 
numerycznej aproksymacji warunków początkowych jest wiec dość 
prosty, zarówno od strony ideowej Jak 1 od strony realizacji 

numerycznej 1 dlatego nie ma potrzeby poświęcać mu wiecej . 
uwagi, 

o wiele bardziej złożony jest natomiast problem 
znalezienia odpowiednich sposobów umożliwiających spełnienie 
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warunków brzegowych w rozwiązaniu numerycznym. Oryginalność 
metody stwarza w tym wzglądzie dodatkowe trudności, gdyż nie 
daje praktycznie możliwości skorzystania z wiedzy 
Literaturowej . 

Przystępując do omówienia tego zagadnienia należy w 
pierwszej kolejności rozstrzygnąć problem podstawowy, dotyczący 
ogólnych możliwości uwzględniania warunków brzegowych w 
rozwiązaniach numerycznych. Istnieją bowiem dwa, różniące się 
dość zasadniczo, sposoby podejścia do tego problemu. 

Pierwszy sposób polega na tworzeniu odrębnych klas punktów 
brzegowych, dla których buduje się specjalne algorytmy 
numeryczne uwzględniające takt, że na brzegu znane są określone 
wartości zmiennych zależnych. Takie postawienie zagadnienia 
komplikuje jednak znacznie cały proces obliczeniowy,. Oprócz 
tego narusza ono jednorodność schematu numerycznego, co często 
prowadzi do generacji niepożądanych numerycznych fluktuacji. Na 
początkowych etapach niniejszej pracy próbowano częściowo 
wprowadzić takie rozwiązanie, a mówiąc konkretnie usiłowano 
przy Jego pomocy spełnić odpowiednie warunki na osi symetrii. 
Wiązało się to z wyodrębnieniem specjalnej klasy punktów 
osiowych, dla których formuły obliczeniowe uwzględniały 
warunki; u= 6Vz-0 dla r~0 oraz zawierały przejścia typu ; 

. Właśnie różnice w numerycznych aproksymacjach 
członu ^ <dla punktów przyosiowych) 1 <dla punktów 

osiowych) były przyczyną trudnych do usunięcia fluktuacji 
parametrów w otoczeniu osi. Dlatego też w późniejszych pracach 
zrezygnowano z wprowadzania punktów osiowych. Oprócz tego 
sposób ten, polegający na bezpośrednim wprowadzaniu wielkości 
wynikających z warunków brzegowych do obliczanych węzłów, 
prowadzi w pewnych przypadkach do trudnych do usunięcia 
numerycznych paradoksów, np, przy modelowaniu procesu 
wychodzenia fali uderzeniowej na swobodną powierzchnię, 

2 wyżej wymienionych powodów zrezygnowano'więc 
definitywnie z tego sposobu i zastąpiono go innym, 
prostszym,bardziej uniwersalnym i nie naruszającym 


73 














Jednorodności schsiuatu obliczeniowego sposobem. Polega on na 
vfprowad 2 eniu tzw, punktów fikcyjnych, którymi nalai:y otoczyć 
brzeg rozważanego obszaru. W punktach tych przechovAijs sie 
wynikające z warunków brzegowych informacje o wartościach 
odpowiednich zmiennych zależnych. Rozważmy wie<^ szczegółowo, 
Jak wykorzystać taką idee <io numerycznej aproksyiuacj i warunków 
brzegowych, nai 

1 . powierzchniach swobodnych, 

2 . powierzchniach swobodnego poślizgu Cgaz - ciało stałe), 

3 . osi syioetril. 

Ad. l. Rozważmy sytuacją pr 2 edstav/lQną na rys. 12. 



Rys.12. Modelowanie warunków brzegowych na powierzchniach 
swobodnych. 


Brzeg badanego ciała tworzą leżące wzdłuż linii S punkty 
A,B,C,D itd. Punkty B,C oraz przykładowo 1,2 1 3 tworzą 
sąsiedztwo punktu A, a punkty B” i C** wynikają z interpolacji 
między parametrami punktów A i B oraz A i C. Punkty 
A’,E*,C',D’, itd, są przyporządkowanymi odpowiednim punktom 
brzegowym punktami fikcyjnymi. 


Jak v/ladonio warunki brzegowe na powierzchni swobodnej 
•sprowadzają sią do zadania na niej sv/lązków: 0 

2 tego co powiedziano powyżej vfynlka, że próby wprowadzania 
tych relacji bezpośrednio do formuł obliczeniowych dla punktów 
brzegowych są nieefektywne. Małeży wiąc zastąpić te relacje 
inną informacją, która byłaby im numerycznie równoważna. Może 
nią być na przykład informacja o stanie ośrodka graniczącego z 
rozważanym obszarem. V pi'zypadku powierzchni swobodne) możemy 
powiedzieć, żo warunek ^*x ® 0 bedzle spełniony na 

powierzchni ciała, Jeśli w ośrodku z nim graniczącym będą 
cerowały sle wszystkie składowe tensora naprążeiir 0 . 

Innymi słowy, warunek w punkcie A możemy zastąpić 

warunkiem 6<ic*0 w punkcie A*. V takiej sytuacji obliczenia w 
punkcie A można prowadzić w oparciu o te same schematy 
numeryczne, które wykorzystuje si^ dla punktów wewnętrznych. 
Zmianie ulegnie jedynie sposób obliczania pochodnych 
przestrzennych, który musi teraz uwzględniać informację zawartą 
w punkcie fikcyjnym A*. 

Praktyczny tok postępowania bedzle wlec następujący? 
a/ V punkcie A wystawiamy normalną do pov/terzchni S 
wykorzystując w tym celu położenia punktów B 1 C, 
b/ V odległości równej promieniowi interpolacji R tworzymy na 
normalnej punkt A‘ i zadajemy w nlmr 6 ^ 4 'jfi. » 0 
c/ Wszystkie pochodne przestrzenne w punkcie A wyliczamy tak 
jakby to był zwykły punkt wewnętrzny, w oparciu o Jego 
punkty sąsiednie: B,C,1,2,3,itd. 
d/ V punktach interpolowanych B** i C” przyjmujemy nowe 
wartości naprężeń, takie, Jakie daje w tych punktach 
płaszczyzna interpolacyjna < 3 . 7 >. 
e/ Kając oki-eślone położenia punktów? A,B'\C“ i A' oraz 

znając w nich składowe tensora naprążert, możemy w punkcie 
A wyliczyć nowe pochodne przestrzenne pól naprężeń. Można 
w tym celu wykorzystać wzory od (3.7) do <3.18) przyjmując 
N=3. Pochodne przestrzenne pola prędkości pozostawiamy bez 
z mia ny, 
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Charakterysxyc2ną cscha taklsgc podejścia jest to, że w 
punkcie A pochodne wzdłuż kierunku normalnego do S są dla 
naprężeń określane na podstawie wartości zadanych na zewnątrz 
obszaru (punkt A'), a dla prędkości na podstawie wartości z 
wnętrza obszaru (punkty 1,2,3,..). Jest rzeczą ciekawą, że 
Inne sposoby aproksymacji pól naprężeń prowadzą do wyników 
obarczonych silnymi fluktuacjami numerycznymi lub nawet o 
niefizycznym charakterze. Bliższe przyjrzenie sią idei ”komórek 
fIkcyJnych"w metodach sieciowych [55,67] wskazuje, źe karzysta| 
si^ tam z analogicznego sposobu obliczania pochodnych 
przestrzennych. Jest on tylko bardziej zamaskowany przez to, 
ciśnienia określa sie tam w środkach, a prędkości na brzegach .| 
komórek. 

Ad.2. Aproksymacja numeryczna warunków na powierzchni 
kontaktu; gazowe produkty detonacji ^ metalowa wkładka 
kumulacyjna okazała sł^ problemem bardziej skomplikowanym. 
Prześledźmy sposób Jego rozwiązania na przykładzie sytuacji 
nokazanej na rys.13. 



o 

O ^ 

produkty 


2 5 n® 

detonac|i 

B 

I / 

° /I 



A 


a ^ 


Rys.13. Modelowanie warunków brzegowych na granicy 
kontaktu: ^ęaz - ciało stałe. 


Powierzchnia S oddziela ciało -stałe od gazav/ych produktów 
detonacji. Punkty A,4,5,6,7,8 reprezentują ciało stałe; B,1,2,3 
-.produkty detonacji; punkt A* jest punktem fikcyjnym 
odpowiadającym punktowi A. 

procedurą zabezpieczającą spełnienie warunków brzegowych dzieli 
sie na każdym kroku czasowym na dwa etapy. Ha każdym z nich 
spełnia slą cząść warunków brzegowych, oddzielni© 1 niezależnie 
dla każdego z ośrodków. 

Veżmy najpierw pod uwagą punkt A reprezentujący brzeg 
ciała stałego. V punkcie tym stosujemy procedurą analogiczną do 
omówionej w poprzednim przypadku. Różnica polega tylko na tym, 
żą w punkcie fikcyjnym A' nie kładziemy ^ i 

p* p ^ 0 . V punkcie fikcyjnym zakodowana jest wiąc 

informa.:Ja, że ciało stał© graniczy z gazem o ciśnieniu p . 
Ciśnienie to na każdym kroku czasowym ustala sią, odpowie-inio 
uśredniając ciśnienia z kilku <np.oznaczonych numerami 1,2 i 3) 
najbliższych wzglądem A' punktów reprezentujących produkty 
detonacji. 

Brzegowy punkt reprezentujący produkty detonacji oznaczono 
na rys.13 literą B. Można przyjąć modelowo, źe gaz 
(reprezentowany przez punkt B) *'wldzi** brzeg ciała stałego 
(reprezentowany przez punkty 4,5,6,7>8 ) jak nieprzenikliwą, 
poruszającą sią ścianką. Rozważmy zatem, w Jaki sposób można 
uwzględnić wpływ tej ścianki na ruch punktu B^. Pole prądkości 
wokół punktu E możemy zamodelować następująco; 

"U = -u o ł- AX-(t- ('f* ) •<* (.-Ł * Zo) (3.22J 

0-= y-e + -O-.j Cł-Zo) <3.23) 

gdzie; , 2. o współ rządne punktu B; - składowe 

wektora prądkości masowej w punkcie B ; ^ Tł / t 

pochodne przestrzenne pola prądkości znajdowane, w oparciu o 
wzory (3.7) (3.18). Pole opisane wzorami (3.22) -5- (3.23) Jest 

polem wokół punktu B pochodzącym od jego punktów sąsiednich w 
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produktach detonacji. Jednocześnie w punkcie B laożna określić 
róv/r 4 leż inne pole prędkości ( 'U- , ) pocłiadzące od 

najbliższych mu punktów brzegowych data stałego <na rys.13 
oznaczonych numerami 4»5,ć,7,8>: 

-u's 14-0 ł- -a'^ ( T-To) + l^-Zo) <3.24 

v’» -iTo + ^,r ~ '*‘ 0 ^ { Z-"Z-o) <3-2E 


i O >Q to fO 

Wielkości: 'U laożna również określić z 

wzorów (3.7) 4 - (3,18) przyjmując, że B Jest punktem centralnym* 
a jego punktami sąsiednimi są punkty: 4,5,6,7,8. Korzystając z 
tych dwóch rozkładów prędkości (. iX ^ '> i < należy 

zbudować trzeci rozkład który uwzględniałby modelowo 

nlsprzenikliwość ścianki oraz możliwość swobodnego ślizgania 
slQ gazu wzdłuż niej. Symbolicznie pokazano te rozkłady na 
rys. 14. 
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Z ”własnego” rozkładu prądkości < tt , 4^* ) tworzymy zatem 
składową styczną nowego pola C 44. f , ) a z 

“zewnętrznego” rozkładu <4t',4r*> składową normalną nowego 

pola. Przyjmując, że w punkcie B brzeg tworzy z osią “z" 


lokalny kąt 

cC , mamy: 


VVt • 

'tLSłVo6 +* \y<jps<C 

<3.26) 

ITk = 

\Ji'coScL ^ 

<3.27) 

Ze składowych i tworzymy nowy, 

prędkości w punkcie B <rys.l4): 

poszukiwany rozkład 

-uf = 

cidióC 4* A^LcosoC 

<3.28) 


CoSoL — 'U"^n.S/4voC 

<3.29) 


Pole ( 44 .^ , ^ ) można przedstawić w postaci; 

'U’' = ^Łe + OtiJ Cr--Tp) + i 4 . 4 Cz.-Z 0 ) 

o-P = m+-('f-fi.) + <r,l(z--2o) 

Po odpowiednich przekształceniach dostajemy wyrażenia na 
pochodne przestrzenne skorygowanego pola prędkości: 

~ -U-itr ł-Co.i^.Cu,'® •** SiW C0S06 < 3 . 32 > 

“^(Z r Co4*oC( 14-4 “' 14 . 12 ) Si«oi coScC ('V^4 * 'Wł) <3-33^ 

■^ir = ('U',1^ - 'U-il) ł- Si-nci CoSoC Cl4.° '-'14-4) <3.34) 

- UJ® 1-Sm'*'c(, (U',1*’*’ U'4) smct CO'ici((l4.,2 - ) (3.35) 


(3.30) 

(3.31) 
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KeasuiDuJekc, w punkcie B pochodne przestrzenne pola prędkości 
oblicza sie ze wzorów <3.32) -t- <3.35)» a pochodne pól napreżeA 
pozostawia sie baz zinlan. V pewnym sensie postępujemy więc 
odwrotnie niż przy obliczeniach w punkcie A, gdzie koryguje słą 
pola naprężeń* a pozostawia bez zmian pola prędkości. 

Tak zaproponcwana procedura daje bardzo dobrą aproksymacją 
warunków brzegowych na powierzchniach kontaktu? gaz - ciało 
stałe. Stwierdzono^ źe procedura ta może okazaó sią zawodną 
tylko wtedy, gdy brzeg ciała stałego nie będzie dostatecznie 
gładki w stosunku do lokalnych gęstości sieci punktów 
brzegowych. Jest to zrozumiałe, gdyż tracą wówczas sens 
operacje u.średnianla pól na wielu punktach, ustalania kierunków 
normalnych i stycznych do brzegu, itp, 

Ad.3» Omawiając niestabilności kątowe zasygnalizowano problem 
konieczności wstawiania punktów fikcyjnych na osi symetrii. 
Wspomniano tam również, że punkty te mogą spełniać podwójną 
rolą, tzn. zapobiegać powstawaniu niestabilności kątowych i 
Jednocześnie zabezpieczać spełnienie odpowiednich warunków na 
osi symetrii. Rozważmy zatem Jakie wartości zmiennych zależnych 
należy przyjmować w osiowych punktach fikcyjnych, takich Jak, 
np.punkt B' łub C* na rys.11. Na osi symetrii spełnione są dwa 
warunki: 'U.-O i 5rz. = 0 i oczywiście te dwie wartości trzeba 
przyjąć w punkcie fikcyjnym. 2 kolei z warunków tych 1 równah 
problemu wynika, że na osi symetrii powinno zachodzić: 



Zerowanie slą tych pochodnych można zasymulować w ton sposób, 
że w punkcie fikcyjnym B* przyjn.jjemy takie same wartości 
funkcji O', p , S<»r 1 Sas. Jak w punkcie B. Reasumując; w punkcie 
fikcyjnym B', którego położenie w każdej chwili czasu wynika ze 
zrzutowania na oś punktu B, należy przyjąć następujące wartości 
zmiennych zależnych: 

<3.37) 
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gdzie indeks "B” oznacza wartości wzlete w punkcie B. Ponieważ 
punkt C <rys.11) nie ma w ogóle przyoslowych punktów sąsiednich 
to odpowiedni punkt fikcyjny C można w metodzie najmniejszych 

iŁwadratów brać z podwójną wagą. 

W praktyce numerycznej często zdarzają, się sytuacje, które 
wymagają dodatkowych zabiegów zabezpieczających rozwiązanie 
przed pojawieniem się nlefIzycznych efektów lub patologicznymi 
deformacjami linii brzegowych. Na rys.15 mamy pokazane dwie 
takie sytuacje wymagające dodatkowych ingerencji. 



Rys,15. Konfiguracje punktów brzegowych wymagające 

dodatkowych działań zabezpieczających stabilność 
obliczeń. 

Rys,16 a)przedstawia sytuację, w której punkt A tworzy 
wierzchołek wąskiego ostrza. Woźna łatwo przewidzieć, że bez 
dodatkowych zabezpieczeń, przy odpowiednio małym cC pojawią 
Sie nlefizyczne fluktuacje parametrów związane z 
niestabilnością płaszczyzny interpolującej pola prądkoścl 1 
naprężeń w kierunku " l "■ Dlatego też wzdłuż tego kierunku 
należy dla dostatecznie małych od dobudować dodatkowe punkty 
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fikcyjne, które na rys. 15 a) oznaczono symbolami A" i A*”. 
Praktycznie, takie zabezpieczenie można wproWadzló już przy 

< 90 '=>. 

Jlatorulast na rys. 15 b) pokazano charakterystyczne 
“wklęśnięcie" linii brzegowej. Takie wklęśnięcia mogą pojawiać: 
się w obszarach, w których zachodzi silne zagęszczanie punktów 
brzegowych. Zagęszczaniu takiemu Eiogą towarzyszyć 
niestabilności wyboczeniowe, ale z reguły ich skala jest mała w 
porównaniu z odległośćiami punktów brzegowych 1 niestabilności 
iizyczne łatwo przeradzają się w niestabilności numeryczne, ' 
których efektem stają się takie, ”jednopunktowe wklęśnięcia". 
Jedynyit skutecznym sposobem przeciwdziałania w takim przypadku 
jest regularyzacja linii brzegowej, a więc usunięcie punktu B z 
obliczeń 1 utworzenie nowego brzegu miedzy punktami 112. 
Działania takie, mające lokalny charakter 1 oczywiście niezbyt 
często podej mowane nie mogą wpłynąć w sposób widoczny na 
globalny wynik symulacji komputerowej, Natomiast ich r 

zaniechanie prowadzi najczęściej do tzw. patologicznych ; 

deformacji linii brzegowych, czyli mówiąc obrazowo do 
"wymieszania się" sąsiednich punktów brzegowych. Oczywiście 
traci wówczas sens procedura stawiania warunków brzegowych i 
cały algorytm ulega rozpadowi. 

3.5. Metoda obliczeń efektów lepko—plastycznych, 

Ogólna postać związków konstytutywnych , z wyodrębnieniem 
członu opisującego Ispko-plastyczne płynięcie materiału może 
być zgodnie z równaniami <2. 15> (2.20> przedstawiona 

następująco; 


Fik. i> SiK 


7si/s,v 


o 


2 ,.. 2 . 


S*/ S./ < -|r Y 


(3.36)=^ 


<3,39) 


(3.40) 
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Jeśli w Jakimś obszarze badanego ciała mamy do czynienia z 
odkształceniami o charakterze sprężystym to zachodzi tam ^=0. 
Jeśli natomiast zachodzi jf>Q to oznacza^ że występują tam 
efekty lepko-^ plastyczne go płynięcia materiału. Zakresy 
odkształceń sprężystych od lepko-plastycznych oddziela warunek 
plastycznego płynięcia Misesat 



<3.41) 


V ciele lepko-plastycznym warunek ten może być dynamicznie 
przekraczany <^>0), a relaksacja naprężeń do stanu opisanego 
warunkiem <3.41) odbywa się w czasie, którego charakterystyczną 
miarą jest, tzw. czas relaksacji X » Przyjmując >^2^1, można 
go w przybliżeniu ocenić; 



Zależność współczynników i p od ciśnienia, gęstości a 
zwłaszcza temperatury powoduje, że w Interesujących nas 
zagadnieniach kumulacji możemy mleć do czynienia z bardzo 
szerokim zakresem zmian czasu relaksacji; HT ~ 10“® + 10"'’* s. 
Tyiaczasera typowe wielkości kroków czasowych Ałiiip* 
zamieszczonych w pracy przykładów) są rzędu 10“® s. Oznacza to, 
że krok czasowy może być w pewnych obszarach znacznie większy 
od-czasu relaksacji. Natomiast Jawny schemat numeryczny typu 
<3. 2) ma sens tylko wtedy gdy spełniony Jest warunek At <<%, 
czyli gdy procesy relaksacji są dostatecznie powolne w 
stosunku do kroku czasowego. Ze względu na czas i efektywność 
obliczeń nie można dążyć do spełnienia tego warunku metodą 
skracania kroku czasowego. Pozostaje więc tylko poszukiwanie. 
innych niż formalny schemat <3.2), sposobów rozwiązania 
problemu. Poniżej przedstawiony zastanie jeden z takich 
sposobów, który uwzględnia fizyczny sens zachodzących zjawisk i 
nie wprowadza numerycznych zaburzeń przy niezmniejszonym kroku At. 
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V plerwszyia etapie obliczamy zmiany naprężeń bez 
uwzględniania efektów lepko-piaatycznych: 

Dla uproszczenia dalszych wzorów przyjmijmy oznaczenie; 

S t n 4 -Ą 

(J odpowiada 

funkcja Y : 


<3.43) 


= 1 - 


■ 5-; 


<3.44) 


V zaleźno-^cl od wartości funkcji możemy mieć do czynienia || 

2 dwoma podstawowymi przypadkami. 

1. Jeśli zachodzi to oznacza, żej , 

a/ pozostajemy w obszarze odkształceń sprężystych Jeśli • .p 
zachodziło -j 

b/ z obszaru odkształceń lepko-plastycznych wchodzimy w 
obszar odkształceń sprężystych nawet bez uwzględniania | 

efektu relaksacji naprężeń jeśli zachodziło ^'^>0. 

Zarówno w wariancie a/ Jak i b/ podstawowy algorytm nie wymaga 

korekty 1 można przyjąć ostatecznie; 

■ 

CL*''*’"* (T ' 

2, Jeśli zachodzi to znaczy, że wchodzimy lub : J 

pozostajemy C^^>0> w obszarze odkształceń lepko-plastycznych. : 
Efekt relaksacji naprężeń można w tyra przypadku ocenić . 

nas tępu J ąco; : - J 


s',“' - s:r‘ ■ 

nt-i 

Odpowiadająca skorygowanym naprężeniom O funkcja 
będzie oczywiście wyrażona wzorem; 


m 


<3.45> - 


m- 


C3.46> 



Jeśli to oznacza, że mimo relaksacji naprężeń 

pozostajemy nadal w obszarze odkształceń lepko-plastycznych i 
możemy przyjąć; Oaic > 

Jeśli natomiast zachodzi to oznacza to, że na kroHu 

czasowym relaksacja naprężeń prowadzi do zejścia poniżej 
granicy plastycznego płynięcia Y- Oczywiście zejście to ma 
charakter formalny, wynikający ze vfzoru <3.45), który nie 
uwzględnia faktu zerowania sle funkcji 0 w momencie osiągania 
granicy plastycznego płynięcia. Mamy wlec w tym przypadku do 
czynienia z sytuacją, w której "naprężenia sprężyste" 
przekraczają granice plastycznego płynięcia <^■''■"''>0), a efekt 
relaksacji sprowadza Je formalnie poniżej tej granicy <fl <0> 
w ramach jednego kroku, czasowego. Oznacza to. iż czas 
relaksacji Jest w tyra przypadku na tyle krótki, że zachowanie 
Sie ośrodka zupełnie dobrze bedzia opisywał model ciała 
spreżysto-plastycznego Prandtla-Reussa. W modelu tyra zakłada 
Sie, że naprężenia mogą osiągaO granice plastycznego płynięcia 
Y,’ale nie mogą JeJ przekraczać., Należy Je wiec skorygować tak. 
aby spełniony był warunek Misesa, Osiąga slą to stosując tzw. 
procedurę sprowadzania do kręgu płynięcia, czyli mnożenia 
wszystkich składowych dewiatora tensora naprężeń przez czynnik 

<1-^’ 


Uł n*' -I 


t ri-f-l 

^5 ł łc 




1 podstawiając ostatecznie; S • Przedstawiona 

powyżej procedura umożliwia w każdym punkcie, w każdej kolejnej 
chwili czasu określanie stanu naprężeń bez wzglądu na relacje 

miedzy czasem relaksacji T . a krokiem czasowym At. 

Odpowiednio do tej procedury należy również zmodyfikować 
sposćb obliczania intensywności deformacji plastycznej. 

Kle można bowiem bezpośrednio korzystać z równania < 2.hi), 
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Jeśli napretżenia określa si^ z modelu Prandtla-Reus^a. Aby 
unikną* tego typu problemów nale±y równanie (3.38) z 
wykorzystaniem równania <2.31> sprowadzić do postaci; 


dt “ Z/4. V di ) 

i dalej wykorzystując (3.43): 




(3.48) 


(3,49> 


Znajomość poszczególnych składowych tensora deformacji 
plastycznej ( & V#J. pozwala bezpośrednio ze wzoru <2.30) 


obliczyć wartość intensywności deformacji plastycznej 


3.6. Modelowanie wzrostu objętości szczelin. 


Veżmy pod uwagę sformułowane w poprzednim rozdziale 
równania opisujące dynamikę powstawania i wzrostu objętości 
szczelin. V równaniach tych występuje progowa wartość średniego 
naprężenia 6 "^« 6 "» , po przekroczeniu której mogą 

zachodzić zmiany objętości szczelin. Okazało się przy tym» że 
wartość 6 'p w pośredni sposób rzutuje na wartość kfoku 
rządowego, przy którym ^^nikl obliczeń nie tracą fizycznego 
sensu. Jeśli przez Ap oznaczymy zmianę ciśnienia na kroku 
czasowym At w jakimś punkcie badanego ciała, to warunkiem 
otrzymania poprawnyvh wyników w jawnym'schemacie numerycznym 
J est; 


A p « e^p 


<3.50) 


Przekształćmy ten warunek tak, aby móc ocenić krok czasowy At 
niezbędny do Jego spełnienia. Z równania stanu (2,21) mamyt 



■I i » 

So 




C3.51> 
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a ze związku (2.36) między gęstością średnią, a gęstością fazy 
ciałostałowej: 


A%s - 8 <,^aVc 


Równanie opisujące dynamikę wzrostu objętości szczelin C2.34) 
pozwala z kolei na uzyskanie następującej oceny <w ramach 
rzędów wielkości!: 


aVc ~ k 6-p Vt At 


Podstawiając <3,53) i (3.52) do (3,51) dostajemy: 


Ap ~ k^gokfff VcA-t 


Warunek <3.50) daje więc następujące oszacowanie. At: 


Uwzględniając konkretne wartości współczynników k|> » k 

oraz zakres typowych zmian 10 iO ^ “ dostajemy: 

Jak widać spełnienie warunku At<< Tc. nie jest możliwe przy 
typowych krokach czasowych A’t ^ 10 ( dla realnych wymiarów 

wkładek kumulacv1 nych pokrytych siecią 10^ -s- 10® punktów). 
Podobnie jak dl.i; efektów lepko-plastycznego płynięcia trzeba i 
w tym przypadku znaleźć takie rozwiązanie problemu, które 
dobrze modelując fizykę zachodzących zjawisk nie pociąga za 
sobą konieczności zmniejszania kroku czasowego. 

Bazując na tym, że w obszarach tworzenia się szczelin 
gęstość fazy ciałostałowej ^5 Jest bliska ^0 , a równanie . 
stanu w tych obszarach ma już przybliżony charakter, można 
przyjąć je w nieco uproszczonej postaci: 


p-kj-l- -|i) -ł-TpSofc 
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Vyicorzyst\jJąc związek < 


P ~ (i-*-f ) + To 806 kn 


<3. 56) 


V schemacie abliczenlowym przed wyliczeniem ciśnienia p, mamy 
Już z równania ciągłości i równania zachowania energii 
wyliczone dla chwili wartości; ^ i . W takim; 

razie wzór <3.58> wiąże na nowym kroku czasowym tylko nieznane: 
wielkości ciśnień i objętości szczelin; 


= Q4 + a*. ( 


gdzie: , Oł - znane na kroku t stałe. V dalszych 

przekształceniach wskaźniki; oraz < u » k. > b^dą 

opuszczane, a wiąc należy pamiętać, że odpowiednie wzory będą 
dotyczyły wybranego punktu dla chwili t Wykorzystując 

wyrażenie C3.59> możemy wyliczyć progowe na kroku ‘t ” * ** 
wartości objętości szczelin: 


a, ł- aj. - 6'a 


V<i, 


a, =-^0 • vt"Vci 


w’"- -(q< ^/(a^•^0t■Vc-^)^-t- 

ZOj 


2La» (3.63) 
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przy czym ; 


'‘O dlq a,7 6'p (3.64: 

Vc, " 0 cilq a^y - 6'f, (3.65: 

Hając obliczone wartości i Vc, można Już prosto określić 

progowe wartości ciśnień na wzrost p lub zmniejszenie 
objętości szczelin; 

Z kolei znajomość wartości progowych i p na zmiany 
objętości szczelin pozwala Już zapostuiować następujący 
algorytm obliczeniowy: 

1. Zakładając, że na kroku czasowym At objętość szczelin nl€ 
uległa zialanle < > obliczamy gęstość fazy 

ciałostałowaj ^ J 


0' ^ 

J 6 “* j \/c^ 

oraz odpowiadające Jej ciśnienie z równania stanu; 


p'""' = ł 0'2*\ Ł"*' ) 


Jeśli spełniany Jest warunek; 


p-<p'-^<p' 


tzn. że znajdujemy się w stanie, w którym ciśnienie Jest niższe 
od wartości progowych na zmiany objętości szczelin. Mażemy więc 
na nowym kroku czasowym przyjąć; 
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2 . 


p = p 


vr^ - 


<3.71) 


Jeśli warunek <3.70) nie Jest spełniony, to musimy 
obliczyć zmianą objątości szczelin na kroku At<lpl>6'p>. 
Vprowadźmy oznaczenie: 


i p/n-i-i 


dla 


<3.72) 

dlo 




Przyjmując za miarą średniego ciśnienia, a wartości p • 
lub za oszacowania ciśnienia progowego Sp w ramacb kroku •; 
czasowego At , możemy równanie < 2.54) przedstawić w postaci : 
różnicowej: 


AVc ^ 


k'p ( Vc + Vco) 


a stąd: 


V (Vt'’ 4- Vc. o) exp {- K y p' ) - Vco 


|| 

■ vii® 


Nowe wartości gęstości fazy ciałostałowej i ciśnienia p ; 

obliczamy analogicznie Jak w punkcie 1: 


p"I 
5 s 


^ rt+^ 


p""' 'f(s'r\ fc”*') 


I 

(3. 70) 


Wartości te możemy przyjąć Jako obliczone n£ aowym kroku 
czasowym: 


P - P ) 


- v^‘ 


z wyjątkiem dwóch przypadków; 
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. j , „ . . ł ^ <+' »^ '4» 

a/ Jeśli zaszło p >p ip <pto należy przyjąć: 




b/ Jeśli natomiast zaszło p l 15”*"“* > p*" , to: 

p =*-p ; Vc “ Vc. <3. 79) 

V przypadkach a/ 1 b/ mamy do czynienia z oscylacjami 
ciśnień wokół wartości progowych. Oscylacje te przy 'Tc. < At 
prowadzą właśnie do -nlefizycznych wyników w jawnym schemacie 
numerycznym. Proponowany model zastępuje te oscylacje petzanien 
wzdłuż wartości progowych <3.7S> i <3,79>. 


3.7. Analiza stabilności schematu różnicowego, dobór kroków 
przestrzennych i czasowych, zewnętrznej dyfuzji 
numerycznej i pseudolepkości, 


Na początku niniejszego rozdziału wspomniano Już, że nie 
ma w chwili obecnej możliwości udowodnienia stabilności 
schematu różnicowego dla nieliniowych układów równań 
różniczkowych cząstkowych stanowiących bazę złożonych symulacji 
komputerowych. Sytuacja taka nie może jednak wpływać hamująca 
na rozwój tego typu prac, z dwóch co najmniej powodów. 

Po pierwsze, potrzeby współczesnych eksperymentów 
wymuszają podejmowanie takich prac bez względu na istnienie lub 
nieistnienie odpowiednich dowodów natury matematycznej . 
Poprawność otrzymywanych rezultatów ocenia się z konieczności w 
sposób ”empiryczny'*, tzn. porównując Je z wynikami 
eksperymentów lub innymi rozwiązaniami teoretycznymi. 

Po drugie, schematy niestabilne eliminują się niejako 
automatycznie z obszarów zainteresowań numeryków - praktyków, 
gdyż generowane przez nie błędy narastają tak silnie i szybko, 
że jakiekolwiek sensowne obliczenia nie mogą być przy ich 
pomocy wykonane. 
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schemat różnicowy typu(3.2) Jest właśnie przykładem takiego 
niestabilnego schematu. Okazało sle Jednakże, że wprowadzając 
niewielkie modyfikacje tego schematu, ograniczone do równań 
zachowania pędu, można uczynić go warunkowo stabilnym. W takiej 
sytuacji celowa staje sie każda, chocl,ażby bardzo uproszczona i 
dotycząca zlinearyzowanej postaci równaA, próba 
przeanalizowania stabilności tego scheinatu. Koźe ona bowiejo 
zasugerować odpowiedź na dwa pytania? 

co należy zrobić dla zapewnienia stabilności schematu 

różnicowego 

czy krok czasowy wynikający z warunku stabilności nie 
bądzle zbyt mały z punktu widzenia technicznej 
efektywności obliczeń ? 

Załóżmy zatem, że ograniczamy sią do analizy płaskiego 
przepływu “cieczy sprężystej" z pominięciem efektów cieplnych! 


S-=-s( 


au . Sir' 


da <3 

S ©'T 

dxr. _ ± ^ 

S 0Ł 

p*o*-(g-5o) 

Podstawiając <3.S3) do (3.801 oraz zakładajac g ^ S® 
dostajemy zlinearyzowaną postać równań wygodną do dalszych 
analiz: 

du _ _ ± ^ c 

dt “ So 

dtv _ i_ 2L <: 

dt §0 


Zgodnie z ogólnym schematem <3.19) można równania (3,84) 
(3.86) zapisać w postaci różnicowej? 



<3.87) 


(3,88) 


<3.89) 


Współczynnik A wprowadzony do równań <3.88) i <3.89) bądzie 
szerzej dyskutowany w dalszych rozważaniach. V tym momencie ’ 
należy tylko zasygnalizować, że symbolizuje on wprowadzoną do 
schematu zewnętrzną dyfuzję numeryczną, warunkującą jego 
stabilność, 

Rozpatrzmy dalej zachowanie się na siatce numerycznej 
przykładowej fourierowskiej składowej zaburzenia: 

-- - 


<3,9o> 

Wykorzystując wzory C3,90) można przykładowo wyliczyć? 

pm£.»pee *■’“ (3.91) 

(3 Q2> 

yuie, 

'Si&co bardziej kłopotliwe jest wyliczenie pochodnych z 
wykorzystaniem wzorów <3.90). Musimy w tym celu przyjąć 
dodatkowe założenie upraszczające, że punkty sąsiednie są 
regularnie rozłożone na okręgu o promieniu ^ A . Wówczas 
pochodne przestrzenne można przedstawić następująca 
(na przykładzie ): 


I""’ 

1 ostatecznie? 
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C3.94> ^ 


Zima cza j iąc ponadto: 


OT T i^,*c. ^Ą*j ^ L-ł^i ^ oCj 


ŁoźeziY równania C3.37> 4 * (3.69) przedstawić w następującej 
postaci: 

ły,*!sin(MA)] o-wi 

^«* 4 ^C(W,i-nioC*)^ U« A - ^ po” Z S*«(itA) <3.97)3 

'a'o'*'''e*^'“^*"’"'^’'^= AT/ A- ^ Po” i siTiCmA) (3.98>, 

Odpowiednia dla tego schematu tzw. macierz przejścia Q jest 


zdefiniowana następująco: 


» GS" 


<3,99>^ 

% 


.^dz i e: 


5 nĄ-A i wł-4 rtł-ł \ 

^[po , U.O , O^o ) 
5" “(po", >Uo" , u-o" ) 

Wprowadzając dodatkov^e oznaczenie: . 

= k.oi^ *■ /moLi. 


(3. 101)'.;^ 


modemy raacierz przejścia na podstawie wzorów (3,90) ■?■ <3.101) J 

przedstawić następująco: 1 

r n 

^ ^Sm(nnA)j 0 j A I 



Warunek stabilności określa sią żądając aby największa wartość 
własna ^ isłacierzy G* była ograniczona następująco: 


-I 


(3.103) 


Odpowiednie równanie na wartości własne macierzy G* ma postać: 

+ .S^^(s,-n'Ł(»i&) + s.nHmA3)l (X~ » 0 < 3 . ao 4 > 

gdzie; 

^ (3.105) 

Warunek (3.103) będzie oczywiście równoważny warunkowi: 


I^I < 1 


<3.106> 


Równanie (3,104) ma trzy pierwiastki; 

= A (3 

; - \A''” ^)^“H •■^‘^'■(sin^(kA) + siii*<mA)’j 

Jeżeli wyrażenie podplerwiastkowe w <3,108) jest dodatnie, 
wówczas warunek (3,106) jest spełniony gdy: 

! X [ X i <3, 

i^W przeciwnym przypadku musi zachodzić: 

< i <3. 




(^in^CKA) + sm^( mA)) ^ i 
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Varunlci <3.109) i <3.111) raożiia ostatecznie zapisać łącznie: 

■*" <3.112) 

Vldać teraz Jasno, że scHeiuat różnicowy typu <3.19) bez 
dodatkowych aiodyf IkacJ 1, czyli z bądzle zawsz- 

Klestabllny. Varunelt <3.112) będzie bowiem niemożliwy do 
spełnienia dla każdego At >0. 

Dla śclsłoScl należy dodać, że w rzeczywistości sytuacja 
Jest nieco bardziej skomplikowana. Schemat <3.19) jest oowlem w 
prezentowanej metodzie nierozerwalnie związany ze sposobem i 
ilością wybieranych punktów sąsiednich oraz aproksymacją 
pochodnych przestrzennych, Vybór 6 punktów sąsiednich < w 
stosunku do 4 dla klasycznej siatki numerycznej) wywołuje już 
określone “uśrednienie" pOl zaburzeń wokół danego punktu czyli 
wprowadza pewną dyfuzję numeryczną, którą można nazwać 
“wewnętrzną". Okazało się jednak, że jest ona zbyt mała dla 
zabezpieczenia stabilności. Z reguły naruszenie stabilności 
rozpoczyna się w obszarach przyoslowych łub obszarach o silnie 
nieregularnym rozkładzie punktów sąsiednich, 

2 tych też powodów należało dodatkowo zmodyfikowali schemat 
numeryczny, wzbogacając go o tzw. "zewnętrzną dyfuzję 
numeryczną". Z warunku <3.112) można wywnioskować,, że należy 
zalngerować w człony i występujące w równaniach 

<3,88) 1 <3.89).Hożna przy tym posłużyć sią następującym 

fizycznym wyobrażeniem o tłumieniu pojawiająch sią fluktuacji,. 
JeSli w punkcie <L,K) wartości zmiennych zależnych zostały 
zaburzone w stosunku do średnich wartości z jego lokalnego 
otoczenia! to można vrprowadzlć automatyczną korektą tego 
zaburzenia zasrąpując np. wartości wyrażeniem: 

Uu.lc Ul.k: ^ A'Cu",k. - U-t.k.) a'- wsp. hcŁbowy <3. 113) 

gdzie: 'Ct oznacza średnią wartość prądkości w punkcie 

;L,K>, wynikającą z prądkości w jego punktach sąsiednich. Tak, 


zaproponowana korekta bądzle prowadzić do '‘wygładzania” pól 
prądkości, a co za tym idzie 1 pozostałych parametrów. 

Takiemu,numerycznemu "wygładzaniu” parametrów, odpowiadałyby w 
układzie równań wyjściowych pewne dodatkowe człony o 
charakterze dyfuzyjnym, i dlatego też efekt takiego 
"wygładzania” nazywa slą cząsto "dyfuzją numeryczną", 

Aby, móc dalej kontynuować analizą warunku stabilności 
C3, 112> należy określić w <3,113) wartość 1 znaleźć 

zależność współczynnika A od A* t A . Zakładając, podobnie 
jak przy określaniu pochodnych przestrzennych, symetrią w 
rozkładzie punktów sąsiednich na promieniu A ♦ można ttL.,K, 
przedstawić nastąpująco: 

r. ^ ^ i <3.114) 

le +e +e -^6 j 

1 dalej po prostych przekształceniach: 


-■w ^ Ai^ ^ ęo$ L»lAl. <3.115) 

Podstawiając (3,115) do <3,113) 1 porównując z <3,38) dostajemy 
wyrażenie dla A : 


\ = ^ (cos(Ra) j- c<ós(mA)) --ł] 


Warunek stabilności <3.112) można z wykorzystaniem zależności 
<3.118) przedstawić nastąpująco: 

X'[i: tcoS.(k^)+cx)s(mA))--t] <3.117) 

X[i:(cos(K/l)+co5(»nA))-l] + £^^(sEnHKA) + S.n^(mA)K0 <3, ua) 

Z warunku <3.117) wynika, że: 


0< i 
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co Jest dość oczywistym zakresem zmian wartości liczbowego 
współczynnika 7k ^ 

Natomiast warunek C3. 118) daje już bezpośrednią możliwość oceny 
kroku czasowego, gwarantującego stabi . :::qść schematu 
numerycznego! 


At < 


asTJP 

2.C 


<3.120) 


Przy 1 warunek C3.120> przechodzi w klasyczny warunek 

stabilności Couranta-Friedrlchsa—Lewy'ego dla schematu Laka. 

V wyniku przeprowadzonych analiz uzyskaliśmy wiąc 
odpowiedz na postawione na wstąpię pytania, odnośnie sposobu 
zabezpieczenia stabilności metody różnicowej oraz wielkości 
kroku czasowego warunkującego tą stabilność. Odpowiedź ta ma 
oczywiście tylko jakościowy charakter i nie wynika z niej 
praktyczny sposób wprowadzania zewnętrznej dyfuzji numerycznej 
na nieregularnej sieci punktów obliczeniowych. Tymczasem Jest 
to problem dość istotny, gdyż dyfuzja ta zabezpieczając 
stabilność 1 tłumiąc numeryczne fluktuacje metody nie może 
równocze -ie zbyt silnie rozmazywać frontów falowych. W ramach 
dotychczas wykonywanych badań najlepsze rezultaty uzyskano 
stosując niżej opisany sposób. 

Załóżmy, że pole prądkóści w lokalnym otoczeniu punktu 
<L,K> opisane Jest następująco: 


u"Kz)= u" ł- + b^Cz-z",jc) < 3 . 121 ) 


\r"(r.z)= u■”^- qa(r-r",»Ł) + ba(z-zC..c) <3.122) 


V odróżnieniu od ekstrapolacji typu <3.7> służącej do 
obliczania pochodnych przestrzennych zakładamy teraz, że 
wartości 1 są nieznane: 
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Ponadto, nie przeprowadzamy interpolacji położeń i parametrów 
na promień R wzglądem punktu <L,K), ale pozostajemy przy ich 
rzeczywistych wartościach. V ten sposób punkt <L,K> nie jest 
wyi‘ćżniony spośród punktów tworzących jego sąsiedztwo. 
Uzasadnienie takiego postępowania wyniknie z dalszych rozważań. 

Utwórzmy zatem wyrażenia; 


Aa.Cir>,,,bO=l 


4. W 

gdzie; 'T* , Zv ** położenia punktów sąsiednich w chwili X , 

wartości składowych prędkości w tych 
punktach, 

]!I a: 9 - w przyjętym systemie wyboru ośmiu punktów 
sąsiednich. 

Współczynniki Cl< , * Ctt » » dL *'* 1 znajdujemy 

przy pomocy metody najmniejszych kwadratów: 

-^=0; 

^--0; 

Wprowadzając następujące oznaczenia: 

o/, -i ^3.129) 

oć^^Ź. -zL.ł) <3.130) 

- 2 ::,^^- <3131) 

^ i =-4 

OĆu = S - /rTiK.) . <3.132) 

C^i 


0 ^5- £ - z.i,k) 

^ i ' 
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z [v" + aj-r."-4- bi(z:" - 2."u.i<)=4 

li-f ^ ^ 

Jeżeli przyjąć, że każdemu z N punktów przyporządkowana Jest 
taka sama masa, to ze związków (3.146> 1 <3.147) wynika, że ped 
punktów rzeczywistych i-U;" . > Jest równy pędowi punktów 

na płaszczyźnie Interpolacyjnej <3.121) i <3.122). Innymi,słowy 
Jeśli każdemu z N punktów dodamy ped; 

ip,"»X[a'’<-a,(ri”-r,';K) *-b,( 2 i«-zS.^)-ui"] <3.>4e> 

to p^d układu N punktów nie ulegnla Zffitanle» bo z <3.146) i 
(3.147) wynika, ż©; 

Współczynnik X nsa ton sam ‘sens co współczynnik wprowadzony do 
związku <3.113). W ten sposób wprowadzona dyfuzja nusueryczna ma 
"nielokalny” charakter, Wprowadzając korekt© prądkości w 
punkcie (L,K) korygujemy równocześnie rozkłady prędkości w 
całym sąsiedztwie tego punktui 


Ui" -H 

> -Ul" + Apt'^ 

( 3 , 151 ) 

\r;'' — 

p- + Zipi"^ 

( 3 . 152 ) 

gdzie! i 

1 nk-Jł lTJVin WVI&affa stOSO’ 

wania relatywnie 


Dyfuzja o charakterze lokalnym wymaga stosowania relatywnie 
większych współczynników 1 niekiedy bywa zawodna, zwłaszcza 

W obszarach przyosiowych. 

pla prezentowanych w pracy przykładoV) 7 ch rozwiązań przyjmowano 
następujące wartości współczynnika A . 
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4 . ĄO-z, 


(3.153)- 


dla materiału vfybucliawQg’0 

A' a 5- lO-*- 


- dla miedzianych wkładek kumulacyjnych, 

Ifa osobna uwagQ zasługuje zaproponowana zależność współczynnika ^ 
od kroku czasowego At, Vynika ona stąd,że w praktyce krok 
czasowy zmienia slq w dość szerokich granicach i należy 
wyeliminować formalne wygładzanie rozkładów prędkości przy At-**0 

(aby zachodziło >u -P- AJ. C,k. £d~*-0 misi 

zachodzić 

Wprowadzona do schematu dyfuzja numeryczna umożliwia 
również modelowanie procesów propagacji frontów fal 
uderzeniowych, rozmywając Je na szerokość kilku kroków 
przestrzennych. Okazało s1q jednakże, że wymaga to stosowania 
większych wartości współczynników A* niż są potrzebne do 
utrzymania stabilności i eliiainacji przypadkowych fluktuacji. 

Aby uniknąć tego niepożądanego efektu, polegającego na zbyt 
silnym "wygładzaniu*' wyników w całym rozważanym obszarze, laożna 
nie podwyższając współczynnika A* , wprowadzić dodatkowo do 
schematu numerycznego pseudolepkość. Wystarcza przy^ tym użycie 
prostej, skalarnej postaci pseudolepkości: 


fcx>nst • 8 (diV Cj )*• dla cl.v CT < O 

dlh. ciiv CiJ 


(3,155) 


Ka zakończenie rozważań zawartych w tym punkcie należy 
Jeszcze omówić problem doboru kroków czasowych i 
przestrzennych. Krok czasowy zapewniający stabilność algorytmu' 
obliczeniowego dobrana na drodze badań testowych, odpowiednio 
rozbudowując Jakościową formułę (3,120)s 


^ O*' 

_3 u,tc ^ S L^ic 
Pu,K. - 


<3,156) 


<3.157) 


I = tS - 0.5 ; 0,1 < I < 1 

gdzie; ~ ciśnienie odniesienia <praktycznle przyjmowano 

P’*= Ph >• 

Punkt (L,K) w którym wyrażenie <3.155) osiąga <po wykonaniu 
każdego cyklu obliczeniowego > najmniejszą wartość 

decyduje o kroku czasowym dla całego algorytmu. 

Krok przestrzenny , rozumiemy Jako średnią odległość 

punktów obliczeniowych, zależy oczywiście od ilości punktów, 
którą możemy pokryć badane ciało uwzględniając moc obliczeniową 
wykorzystywanej EMC. V dotychczasowych badaniach z 
wykorzystaniem EMC R-60, liczba punktów obliczeniowych nie 
przekraczała 1500 (łącznie dla materiału wybuchowego 1 
wkładki). Odległość najmniejszego zbliżenia ^ mfn przyjmowano 
bliską 0,25 R. . 

Celem zminlmaltzovfania ilości ingerencji w proces 
obliczeniowy, wynikających z eliminacji niestabilności 
kątowych, stosowano sieć przestrzenną specjalnie zagęszczaną 
wzdłuż linii brzegowych. Efekt ten jest dobrze widoczny na 
rysunkach zamieszczonych w następnym rozdziale, a 
przedstawiających powiększone sieci numeryczne, 

3.8. Organizacja procesu obliczeniowego. 

V ramach dotychczasowych rozważań zawartych w niniejszym 
rozdziale, przedstawiono koncepcję numerycznego rozwiązania 
układu równań różniczkowych opisujących sprzężone procesy 
detonacji materiałów wybuchowych oraz deformowania 1 napędzania 
przez nie ciał stałych. Aby móc tę koncepcję zweryfikować, 
ocenić JeJ wady i zalety, zbadać zakresy zastosowań metody, a w 
dalszej kolejności uzyskiwać rozwiązania konkretnych problemów, 
należało zbudować odpowiedni dla tej koncepcji kod komputerowy. 
Nazwą "kod komputerowy" przyjęto określać w literaturze' 
programy komputerowe służące do rozwiązywania wielowymiarowych, 
niestacjonarnych problemów fizyki matematycznej , ze względu 








ua ich złożoność od strony matematyczno-fizycznej oraz 
numerycznej i technicznej. 

Właśnie problemom natury technicznej, a wlec związanym z 
organizacją procesu obliczeniowego 1 budową samego programu, 
należałoby poświecić Jeszcze nieco uwagi, Niestandardowośó 
metody pociągnęła bowiem za sobą konieczność stosowania 
specyficznych rozwiązań technicznych. Rozwiązania te omówione 
zostaną ogólnie i tylko od strony ideowej., gdyż sam program 1 
Jego eksploatacja to odrębny, dość złażony i Jednocześnie zbyt 
szczegółowy problem, aby włączać go w ramy rozważań niniejszej 
pracy. 

Podstawowe znaczenie dla funkcjonowania kodu numerycznego 
ma problem przechowywania, przesyłania oraz uaktualniania 
Informacji o stanie ośrodka w każdym punkcie obliczeniowym, 

ITaJprościej problem ten można rozwiązać umieszczając te 
informacje w odpowiednich tablicach, dla których rezerwuje sle 
określone obszary pamięci operacyjnej komputera. Tablice te, ze 
względu na ich funkcje, można podzielić na cztery grupy. 

Do pierwszej grupy należy zaliczyć tablice służące do 
przechowywania wartości zmiennych zależnych w każdym punkcie, w 
chwili t i t . Kają one postaćr PTS(L,K,N> i PT<L,K,N). 

Punkty obliczeniowe pokrywające badany obszar ponumerowano 
dzieląc Je na linie 1 porządkując je wzdłuż tych linii. Para 
liczb CL,K> oznacza więc K-ty punkt z L-teJ linii i może być 
traktowana jako Jego lagrange*owskle współrzędne. V trzeciej 
kolumnie przechowywane są wartości zmiennych zależnych, np. 

PTSCL.K.i) = •, PT(U,KW)=r:rK 0 

PTS (L, (C.S) = ^ PT (L,K,5> (3.158) 

PTS (L,K,5)= , 

Wymiary tych tablic są zdeterminowane liczbą punktów 
obliczeń. ^ch i ilością zmiennych zależnych. Po wykonaniu 
obliczeń na kroku następuje podstawienie wartości z 


tablicy PT do tablicy PTS, a do tablicy PT można Już wpisać 
nowe wartości z kroku t . 

Do drugiej grupy należy zaliczyć dwie tablice, w których 
przechowywane są informacje □ numerach punktów tworzących 
sąsiedztwo punktu obliczanego. Pierwsza z nich ma 
postać:NDKTCL,K,N), a Informację w niej zakodowaną należy 
rozumieć następująco: 

NDKT(L,K,‘1)=L1 NDK.T(L,K,5)'L2. 

<3.159) 

NDKT NDKT(L,K.,4)=K2. 

Punkty o numerach <L1,K1), <L2,K2>, Itd,tworzą najbliższe 

sąsiedztwo punktu <L,K>. Dzięki tej tablicy każdy punkt 
"pamięta" numery swych punktów sąsiednich, a więc i 
automatycznie Ich parametry, np,położenie punktu CL1,K1> 
znajdujemy następująco: 

r, " « PTS [ NDKT (L.K .i ) , NDKT (L,K,2.), 1 ] 

<3.160> 

Z./^PTSLnDK.tCU.K.-I), NDKT(L,K,2.J, 2.1 

Proces wymiany punktów sąsiednich zostaje w kodzie numerycznym 
zarejestrowany przez zmianę numerów punktów pamiętanych w 
trzeciej kolumnie tablicy NDKT. 

Problem numerycznej aproksymacji warunków brzegowych 
wymaga aby punkty tworzące brzeg były odróźnialne od punktów 
wewnętrznych, a ponadto aby każdy punkt brzegowy "pamiętał" 
najbliższe .sobie punkty brzegowe. Vymaganla te można spełnić 
wprowadzając tablicę "sąsiednich punktów brzegowych" typu: 
NBRZ<L,K,N>. W trzeciej kolumnie tej tablicy zapamiętywane są 
numery identyfikacyjne "sąsiadów brzegowych" analogicznie Jak w 
relacjach <3.159), Punkty wewnętrzne można dzięki tej tablicy 
łatwo odróżniać od brzegowych,gdyż dla punktu brzegowego 
NBRZ<L, K, 1 >=L1?^0, a gdy punkt <L,K) Jest punktem wewnętrznym, 
to można np. zadać NBKZ(L,K,1>“0, Zawartość tablicy NBRZ Jest 
korygowana, gdy zachodzi wyeliminowanie z obliczeń punktu 
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brzegowego, z-aiulana punktu v/evmQtrzneso na brzegowy lub przy 
modelowaniu fragmentowania slsj ośrodka. 

Trzecią grupą tablic tworzą tablice pomocnicze, w których 
zapamiętywane są inforroacje nie mające tak podstawowego 
znaczenia dla procesu obliczeniowego, ale przydatne z różnych 
technicznych powodów. Przykładowo, są to Informacje o; 

-czasie dojścia frontu fali detonaoyjnej do poszczególnych 
punktów reprezentujących materiał wybuchowy, 

-położeniach punktów fikcyjnych związanych z 
poszczególnymi punktami brzegowymi, 

-poprawkach do rozkładów prędkości wnoszonych przez 
zeymątrzną dyfuzją numeryczną. 

Wreszcie, czwarta grupa tablic, to tablice o wyraźnie 
pomocniczym charakterze, służące do chwilowego prżechQwyv/anla 
pewnych wielkośćl, aby uniknąó ich wielokrotnego 
obliczania itp. 

Kod komputerowy został napisany w jązyku FORTRAN i ma 
typową dla tego jązyka strukturą - składa slą z programu • 
głównego 1 wywoływanych z niego podprogramów. Informacje o 
zawartości wymienionych wyżąj tablic oraz o wartościach różnych 
współczynników są przekazywane do podprogramów za pomocą 
instrukcji COMJiON. 

Scharakteryzujmy zatem krótko podstawowe funkcje jakie 
spełnia program główny i poszczególne, vfażniejsze podprogramy. 
Symbole tablic i podprogramów odpowiadają symbolom używanym w 
kodzie numerycznym, 

1. Program główny CHAIN). 

a/ Deklaracja struktur podstawowych tablic, 

b/ Wprowadzenie cząścl danych początkowych i współczynników 
określających: 

- wielkości kroków przestrzennych 1 czasowych, 

“ częstości 1 sposoby wywoływania podprogramów, 

- wymiary sieci numerycznych dla materiału wybuchowego i 
wkładki kumulacyjnej, 

- róvmania stanu,model Steinberga, wartości w punkcie C-J, 


lOó 


c/ 

d/ 

0 / 

f/ 

g/ 

2 . 

a/ 

b/ 

c/ 

d/ 

0 / 

3. 

a/ 

b/ 

c/ 

d/ 


wielkości zewnątrznej dyfuzji numerycznej i 
pseudolepkośći, 

“ sposoby sterowania pracą całego programu. 

Wywoływanie poszczególnych podprogramów. 

Wymiana zawartości tablic typu: PTS=PT, po wykonaniu 
pełnego cyklu obliczeniowego na danym kroku czasowym. 

Dobór kroku czasowego. 

Korekta rozkładów prądkośoi związana z zewnątrzą dyfuzją 
numeryczną. 

Organizacja wydruków 1 zapisów wyników pośrednich w 
parniącl dyskowej. 

Podprogram wyboru punktów sąsiednich <SAS>. 

Ustalenie zbiorów punktów sąsiednich w chwili początkowej. 
Tworzenie dla każdego punktu podzbioru najbliższych 
punktów sąsiednich. 

Eliminacja z podzbioru najbliższych punktów sąsiednich 
punktów nie mogących tworzyó fizycznego sąsiedztwa z danym 
punktem (podrozdział 3.2). 

Podział podzbioru najbliższych punktów sąsiednich na 
sektory kątowe. 

Wybór najbliższych punktów sąsiednich z poszczególnych 
sektorów. 

Podprogram obliczania wartości zmiennych zależnych po 
upływie kroku czasowego CPTZW). 

Ustalenie położeń 1 wartości zmiennych zależnych w 
punktach sąsiednich na płaszczyźnie czasowej t , 
Interpolacja położeń i parametrów punktów sąsiednich na 
promień R. 

Obliczanie pochodnych przestrzennych składowych pół 
prędkości i naprężeń. 

Ustalanie położeń punktów fikcyjnych dla punktów 
brzegowych 1 korekta pochodnych przestrzennych wynikająca 
z wartości brzegowych zadawanych w punktach fikcyjnych. 
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e/ Obliczanie pochadnycJi 1 śFsdnich wartości pól prędkości 
związane z zewnętrzną dyfuzją numeryczną, 
f/ Wyszukiwanie najbliższych punktów sąsiednich z drugiego 
ośrodka na powierzchni kontaktu; produkty detonacji - 
wkładka kumulacyjna. Konstrukcja wzajemnych pól prędkości 
i ciśnie* dla realizacji warunku brzegowego na tej 
powierzchni. 

g/ Obliczanie efektów związanych z lepko-plastycznym 
płynięciem materiału i tworzeniem sie szczelin w/g 
specjalnych (podrozdział 3,5 i 3,6>. algorytmów. 

h/ Ostateczne wyliczanie nowych wartości zmiennych zależnych - 

t »n- -ł 

■ 

4. Podprogram służący do budowy początkowej sieci punktów 
obliczeniowych i zadania warunków początkowych (VARP1>. 
a/ Budowa sieci położę* początkowych dla punktów 
reprezentujących materiał wybuchowy i wkładką, 
b/ Wprowadzenie wartości odpowiadających warunkom początkowym 
do odpowiednich- tablic. 

c/ Konstrukcja tablicy "sąsiedztwa brzegowego" - HBRZ. 

5. Podprogram elIrainacJi niestabilności zbliżeniowych 
(BRZEG). 

a/ Eliminacja punktu wewnętrznego przy niestabilności 

zbliżeniowej typu: punkt wewnętrzny - punkt wewnętrzny, 
b/ Eliminacja punktu v/ewnetrznegQ przy niestabilności 

zbliżeniowej typu; punkt wewnętrzny - punkt brzegowy, 
c/ Eliminacja punktu brzegowego przy niestabilności 

zbliżeniowej typu: punkt brzegowy - punkt brzegowy 
(korekta tablicy NBRZ>. 

6, Podprogram eliminacji niestabilności kątowych (BRZEGI), 
a/ Kontrola sieci numerycznej dla identyfikacji punktów 
zagrożonych możliwością pojawienia sie niestabilności 
kątowych. 

b/ Zamiana punktów wewnętrznych na punkty brzegowe dla 

uniknięcia niestabilności kątowych (korekta tablicy NBRZ), 


lOS 


•7. 


Pcciprogram eleminacjl punktów brzegowych w obszarach 
zągrożonych patologicznymi deformacjami kształtów linii 
brzegowych (BRZĘK). 
ą/ Kontrola stanu deformacji linii brzegowyc-h, 
b/ Usuwanie z obliczeń punktów zagrożonych nlestabllnościami 
(korekta tablicy NBRZ). 

'|<llmo» że wymieniono tylko (abstrahując od sposbów 
.realizacji) najważniejsze podprogramy i ich podstawowe funkcje, 
widać, iż kod komputerowy różni sie znacznie od "programu 
komputerowego" w potocznym tego słowa znaczeniu. Kod 
komputerowy Jest zespołem wielu, czysto dosyć złożonych 
pro.gramów wzajemnie ze sobą współpracujących i realizujących 
poszczególne elementy składowe pełnego modelu numerycznego 
badanego zj awiska. 

Tymi uwagami można już zakończyć omawianie kompleksu 
problemów związanych z matematyczno — fizyczno “ numerycznym 
modelowaniem zjawisk kumulacyjnych i przejść do omówienia 
przykładowych wyników symulacji komputerowych. 









Rozdział 4 


PRZYKŁADOWE WYNIKI SYHULACJ I K014PUTER0VYCH 

Przedstawione w poprzednich rozdziałach równania oraz 
sposób ich numerycznego rozwiązania mogą byó wykorzystane do 
modelowania komputerowego bardzo szerokiej klasy zjawisk 
fizycznych dających sie opisać równaniami mechaniki ośrodków 
ciągłych, Cel niniejszej pracy spcyjodował oczywiście 
ograniczenie .slą w dotychczasowych badaniach do klasy zjawisk 
kumulacyjnych 1 tych też zjawisk bądą dotyczyć głównie 
prezentowane rozwiązania. Jednakże, aby udowodnić szersze 
możliwości zastosowań proponowanej metody, zademonstrowano w 
trzecim punkcie niniejszego rozdziału przykład dotyczący 
modelowania przebijania osłony ciałem odkształcaInym. Wybór 
tego zjawiska był nieprzypadkowy, gdyż zagadnienia kumulacji 1 
przebijania osłon muszą być w przyszłości traktowane łącznie, 
aby móc w pełni ocenić efektywność badanych układów 
kumul^cyj nych, 

V rozdziale niniejszym wyniki obliczeń teoretycznych' 
zostaną skonfrontowane z wynikami specjalnie w tym celu 
wykonanych w IFPILM eksperymentów. Dokonane zostaną również 
odpowiednie porównania wyników obliczeń numerycznych z wynikami 
które można otrzymać w oparciu o klasyczną teorie kumulacji lub 
wynikami teoretyczno - eksperymentalnymi zawartymi w 
11teraturz^. 

4.1. "Kumulacja ody/rotna", 


zjawisko "kumulacji odwrotnej" (wybuchowego kształtowania 
pociskowi costało od strony jakościowej omówione w rozdziale I 
Vspomniane tam również, że zjawisko to zainspirowało podjecie 
tematu niniejszej pracy 1 stanowiło podstawę, na której 
budowano zarówno model matematyczno - fizyczny Jak 1 metode 
oraz kod numeryczny 13.261. Weryfikacje wyników teoretycznych 1 


eksperyiuentalnych prowadzono w oparciu o iotografie wybuchowo 
ukształtowanych pocisków, których charakterystyki (kształty, 
wymiary, prędkości) zależały od aktualnych celów i potrzeb 
niezależnie prowadzonych eksperymentów. Natomiast przykład 
przedstawiony w pracy niniejszej dotyczy układu specjalnie 
zaprojektowanego (wymiary, masa materiału wybuchowego) z myślą 
0 teor.etyczno - eksperymentalnej weryfikacji wyników. Jego 
Wymiary pozwalają bowiem, aby oprócz rejestracji optycznych, 
można było w przyszłości wykonać dla niego odpowiednie badania 
Z wykorzystaniem laboratoryjnych urządzeń rentgenowskich. 

Kształt oraz wymiary tego modelowego układu do badania 
żjawlska “kumulacji odwrotnej" przedstawiono na pierwszym 
kadrze rysunku nr 16. 

średnica wkładki wynosi 3 cm, grubość 0.26 cm, a kąt rozwarcia 
2 =150®, Na tym oraz na dalszych rysunkach tego typu 
naniesiono podziałką, której jednostka odpowiada zawsze 
długości równej - 1 cm. Do konstrukcji układu wykorzystano 
termoplastyczny loateriał wybuchov/y, którego właściwości fizyko¬ 
chemiczne opisano równaniem stanu typu (2,7) 1 przyjęto; 

CM.8- ^ ; cT* 0.5 

Punkty, z których utworzone są poszczególne kadry na 
rys.16, są wykorzystywanymi w metodzie numerycznej materialnymi 
punktami ośrodka, tworzącymi odpowiednią sieć numeryczną. 
Zgodnie z przyjętą metodą przemieszczają sie oae wraz z 
ośrodkiem, a wlec dobrze Ilustrują proces jego deformowania 
ślą. Kolejne kadry na rys,16 zostały wykonane w stałych 
odstępach czasowych <2 ps>. Na pierwszym kadrze <t=0> 
zaznaczono strzałką punkt, który przyjęto Jako centrum 
inicjacji sferycznej fali detonacyjnej. Ostatni kadr ilustruje 
sytuacje, w momencie, w którym ciśnienie produktów detonaoj w 
żadnym punkcie nie przekracza wartości około 5 kbar. 

Z licznych eksperymentów komputerowych wynika, że w takim 
przypadku można już w dalszych obliczeniach nie uwzględniać 
wpływu produktów detonacji na wkładkę kumulacyjną. 
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Wy 0 li minowanie z dalszych obliczeń produktów detonacji ina duże 
znaczenie praktyczne, gdyż umożliwia znaczne skrócenie t 
uproszczenie tych obliczeń. 

Ilustracja przedstawione na rys.10 miały na celu pokazanie 
przede wszystkim <ze wzglądu na skale przestrzenne zjawisk) 
kształtów rozprąi-ającego sią obłoku produktów detonacji. Aby 
dokładnie prześledzić zmiany kształtu wkładki oraz zilustrować 
dalsze fazy zjawiska pokazano na rys.17, 18 i 19 deformacje 
sieci punktów w dużej skali przestrzennej 1 z większymi 
odstępami czasowymi miądzy kadrami (10 ps) . 

Ifa rys. 17 widzimy wkładką w stanie niezaburzonym (pierwszy 
kadr) oraz po napądzenlu 1 wstąpnej deformacji (drugi kadr). 
Drugi kadr z rys.17 oraz pierwszy kadr z rys.16 pokazuJa> ż® 
badany proces ma w tej fazie pewne cechy "odwrotnego procesu 
kumlacji". Wkładka jest bowiem nachylona do osi pod kątem, 
który w przeciwieństwie do kumulacji klasycznej, jest wiąkszy 
od 90®. Ponadto wystąpują tam elen^nty przypominające zbitek 1 
strumień, z'Jednak różnicą*, że w kierunku przegrody (z lewa 
na prawo) pierwszy podąża zbitek, a dopiero za nim szczątkowy 
strumień. Okazuje sią jednakże, że dla dalszych chwil czasu 
analogie z klasyczną teorią kumulacji już nie wystąpują. 
Hastąpuje bovd.em wyrównanie prądkości "zbitka" 1 "strumienia" i 
zakończenie procesu deformowania sią wkładki w otoczeniu osi. 

Z pewnym opóźnieniem zostaje wyhamowany ruch ramion wkładki w 
kierunku osi i wyrównanie prądkoścl poosiowej, V rezultacie, w 
kierunku przegrody porusza sią jedno ciało, którego kształt nie 
zmienia sią Już w dalszych chwilach czasu <U=0, V=CQnst). 
Kształt ten widoczny jest na drugim kadrze rysunku 19. 

Okazuje sią wiąc, że omawiane zjawisko ma jakby dwie fazy: 
hydrodynamiczną 1 wytrzymałościową. V fazie hydrodynamicznej 
wystąpują stosunkowo duże gradienty prądkości i znaczne siły 
bezwładności w stosunku do sił wytrzymałościowych 
< 5» V ). Wkładka kumulacyjna ma w tej fazie własności 
podobne do cieczy i to jest powodem, że zachowuje sią ona w 
sposób zgodny Jakościowo z klasyczną teorią kumulacji. 




































Gradienty pr^dko^cl przy wkładkach o dużych kątach rozwarcia 
nie są Jedak na tyle duże, aby po pewnym czasie nie zostały 
wyrównane pod wpływem sił wytrzymałościowych. Faza 
hydrodynamiczna przechodzi wtąc w sposób ciągły w fazą 
wytrzymałościową, w której efekty sprężysto lepko - 
plastyczne decydują o ostatecznycznym kształcie zdeformowanej 
wkładki. Rozwiązanie komputerowe dostarcza wlec dodatkowego 
argumentu na rzecz zastępowania w przyszłości terminu "odwrotna 
kumulacja" odpowiedniego tylko dla hydrodynamicznej fazy 
zjawiska, terminem “wybuchowe kształtowanie 1 napędzanie 
pocisków" obejmującym obie fazy procesu. 

Analizując poszczególne kadry przedstawione na rys.17+19 
można sle zorientować, że dwukrotnie w trakćie obliczeń 
przeprowadzono operacj e oderwania brzegowych fragmentów 
wkładki. Należy'wląc udzielić odpowiedzi na kilka pytań 
dotyczących zaróvmo celów jak 1 sposobów przeprowadzania takich 
operacji. 

V trakcie obliczeń okazało sie, że dla układów o geometrii 
ładunku 1 wkładki zbliżonej do przedstawianych na rys.16 
występuję zawsze efekt odrywania sią brzegów wkładki od jej 
zasadniczej części. Kożna przy tym wyróżnić Jakby dwie fazy 
tego zjawiska. V pierwszej fazie stosunkowo niewielka msa 
tworzącą brzeg wkładki uzyskuje pod wpływem ciśnienia produktów 
detonacji prędkość U >0,a wiec ekspanduje na zewnątrz układu. 
Prędkość tej ekspansji jest zbyt duża aby siły wytrzymałościowe 
były zdolne zapobiec oderwaniu sle tego fragmentu. Katomiast w 
drugiej fazie ma miejsce bardziej długotrwały proces dalszego 
fragmentowania sie brzegów wkładki, spowodowany zbyt małą 
prędkością poosiową tych fragmentów w stosunku do prędkości 
zasadniczej części wkładki. 

Operacja wyeliminowania z obliczeń fragmentu, który 
fizycznie oderwał slą od wkładki. Jest niezbędna Jeśli 
oczekujemy zgodności wyników teoretycznych 1 eksperymentalnych 
dla dalszych faz zjawiska. Wynika to stąd, ie Jeśli fragment 
taki nie zastanie wyeliminowany z obliczeń, to pochodne 
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<4, 3> 


przestrzenna bądą liczone dalej w oparciu o punkty 
reprezentujące ten fragment, a wiec punkty wkładki sąsiadujące 
z tym fragmentem będą się znajdowały w niefizycznych polach 
naprtjźeń 1 pr^dko^ści, Rozbleźno^ści między eksperymentem 
komputerowym a eksperymentem fizycznym będą w takiej sytuacji 
narastały nawet v/6wczas, gdy fragment taki mo±na uznać za obłok 
nieoddziaływujących odłamków <6'^’ic.= 0, 

Załóżmy chwilowo, że potrafimy wskazać obszary (punkty), w 
których nastąpiło zniszczenie materiału. Operacje oddzielenia 
fragmentu od wkładki można przeprowadzić tylko wówczas, gdy 
zniszczenie materiału nastąpiło w obszarze całego przekroju 
poprzecznego ramion wkładki. Punkty reprezentując© oderwany 
fragment eliminuje się z dalszych obliczeń, Dla pozostałych 
punktów ustala się nowe tablice sąsiedztwa 1 sąsiedztv/a 
brzegowego i obliczenia mogą być dalej kontynuowane, 

Do omówienia, w ramach tego problemu, pozostało Jeszcze 
zagadnienie kryteriów Jakie były przyjiaov/ane dla określenia 
miejsc, w których nastąpiło zniszczenie materiału, PQnlQv/aż w 
zależności od czynników zewnętrznych oraz v/łasnoócl 
materlało^wych zniszczenie może mieć charakter ctągliwy (duża 
deformacja plastyczna v/evmątrz ziaren), kruchy (wzrost 

objętości szczelin w przestrzeniach między ziarnami) lub 
najczęściej mieszany, wprowadzono dwa kryteria zniszczenia. 

Jedno iest określone przez krytyczną wartość deformacji 
plastycznej 6. , a drugie przez krytyczną wartość objętości 

szczelin (Vc,) . Jeśli w jakimś punkcie ośrodka zaszło; 

lub Vc>(Vc)^^ <4.2> 

to przyjmowano, że w tym punkcie nastąpiło zniszczenie 
® *" i a ł u, JJa podstawie badań własnych oraz sugestii 
literaturowych przyjmowano dla miedzianych wkładek 
kumulacyjnych w obszarach niezbyt wysokich temperatur (typowe 
temperatury w obszarach fragmentacjl są rzędu (500 e 600 K>: 



Ka rys,20 przedstawiono dla omawianego przykładu sytuację 
w momencie pierwszej, a na rys,21 drugiej fragmentacji, Ha obu 
rysunkach zaznaczono linię, którą przyjęto za oddzielającą 
oderwany fragmnt od pozostałej części wkładki oraz naniesiono w- 
jej otoczeniu wartości 1 *Ve. . Widać, źe w otoczeniu linii 

fragmentacjl zachodzi na rys.20: 

£,r > t J,. ; powierzchni wewnętrznej (od 

strony materiału wybuchowego), 

ŁP ^ i "VŁ>(Ve) - na powierzchni zewnętrznej , 

a na rys.21: 

tp > t. » M<(Vc)k, - na powierzchni wewnętrznej, 

ep).ep^ i Vc>(X)icf- na powierzchni ze\*mętrznej , 

Z analizy tych relacji oraz relacji geometrycznych 
uwidocznionych na rys.20 1 21 wynika, że jeśli kryteria 
zniszczenia materiału mają być spełnione z pewnym "zapasem”, to 
pozostaje określony margines swobody odnośnie wyboru linii 
fragmentacjl oraz czasu przeprowadzenia tej operacji. 
Eksperymenty komputerowe pokazały, że margines tan Jest dość 
szeroki, tzn, nie obserwuje się istotnych zmian parametrów 
formowanego pocisku Jeśli linia fragmentacjl będzie 
poprowadzona przez inne punkty niż to zaznaczono na rys,20 i 21 
z dokładnością do 3 punktów, Czas przeprowadzenia operacji 
rozdzielenia rozwiązań na dwa obszary może być również znacznie 
opóźniony w stosunku do przyjętego w tym przykładzie, nawet w 
zakresie do 10 ps, 

Szczegółowa obserwacja zmian wielkości i *Vc w 

obszarach fragmentacjl wskazuje, że kryteria zniszczenia 
materiału są spełniane stopniowo dla kolejnych punktów 
począwszy od podstawy wkładki. Oznacza to, że fizycznie mamy do 
czynienia z procesem ąuaslciągłego odrywania się odłamków przy 
podstawie wkładki. numeryczne prowadzenie takiego ąuaslciągłego 
procesu fragmentacjl Jest formalnie możliwe, ale technicznie 
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Rys.20.Rozkłady deformacji plastycznej oraz objętości 

szczelin Vc vf chwili przeprowadzania pierwszej 
operacji, oddzielenia zniszczonego fragmentu 
wkładki, 

bardzo kłopotliwe. Dlatego też bazując na tym, że czas 
przeprowadzania operacji rozdzielania rozwiązań może być 
przyjmowany z dożd dużym opóźnieniem, podzielono proces 
fragmentacji tylko na dwa, wyżej omówione etapy. 

Fragmentacja obszarów przy podstawie wkładki kumulacyjnej 
j^st efektem potwierdzonym eksperymentalnie. Je^ll odległość 
ładunku kumulacyjnego od przegrody Jest niewielka, to można 
obserwować radialnie rozłożone ślady po uderzeniach tych 
fragmentów w przegrodą. Niekiedy można również zaobserwov/ać te 
odłamki na fotografiach wybuchowo kształtowanych pocisków. Na 
rys.22 zamieszczono szczególnie ciekawą fotografią. Dobór 
parametrów ładunku spowodował bowiem oderwanie sią od wkładki 
dużego fragmentu, który nie rozproszył slą radialnie lecz z 
nieco mniejszą prądkością podąża za pociskiem. 
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Rys. 21.Rozkłady deformacji plastycznej oraz 

objątości szczelin Yć. w chwili przeprowadzania 
drugiej operacji .oddzielenia zniszczonego 
fragmentu wkładki. 


Zgodnie z tym co powiedziano we wstąpię do niniejszego 
rozdziału, rozwiązanie komputerowe pozwala na sporządzenie 
bardzo wielu różnych czasowo - przestrzennych charakterystyk 
procesu. 

Na rys.23 przedstawiono kilka takich przykładowych 
charakterystyk. Dotyczą one rozkładów temperatur, średnich 
gąstości 1 deformacji plastycznych wzdłuż odpowiednich 
przekrojów uformowanego ostatecznie pocisku oraz zmian w czasie 
prądkości czoła ( ^ ^ tylnej cząści pocisku ^ • 

Rozkłady gąsxoścl wzdłuż osi <przekrój A—A’), oraz w 
zasadniczej, czołowej cząści pocisku CprzekróJ B-B') wskazują, 
że w tych obszarach ośrodek ma strukturą litego metalu z 
pomijalnle małą objątością szczelin. Wyjątek stanowi tu tylko 
niewielkie otoczenie punktu A*, gdzie bliski jest spełnienia 
warunek pojawia sią dopiero w końcowych 
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Rys.22> Fotografia Y^bucbowo ukształtowanego pocisku wraz 
z oderwanym od niego fragmentem podstawy wkładki, 

fazach procesu wyhamowania pr^dko^ól czoła pocisku 1)^ » tj. w 

czasach rzqdu 25 30 ps i Jest on wynikiem tego, że w procesie 

napędzania i deformacji wkładki otoczenie punktu A' osiągnęło 
prędkości bliskie krytycznym. Gdyby konstrukcja ładunku 
umożliwiła osiągnięcie większych prędkoScl 1/ 
fragment czoła mOgłby oderwać się od zasadniczej części 
pocisku, gdyż siły wytrzymałościowe byłyby Już zbyt małe dla 
zrównoważenia sił bezwładności. Efekt taki zarejestrowano 
wielokrotnie na fotografiach wybuchowo kształtowanych pocisków. 
Przykład takiej fotografii zamieszczono na rys.24. 

Efekt odrywania się czołov/ych fragmentów pocisków 
(Yc. >(Vc) jtT pewnym otoczeniu czoła pocisku) występował 

również przy analizach numerycznych różnych wariantów układów 
kumulacyjnych. Dlatego też kod numeryczny Jest przygotowany na 
przeprowadzenie operacji oderwania i wyeliminowania z dalszych 
obliczeń takiego fragmentu. 


122 





















Rys.24. Fotografia wybuchowo ukształtowanego pocisku wraz 
2 oderwanym fragmentem z jego czoła. 

Ideowo zasady postępowania sa w takim przypadku analogiczne Jak 
dla fragmentujących sie ramion wkładki. 

Rozkład temperat»jry wzdłuż osi (A-A‘ > przedstawiono 
dlatego, że na osi występują temperatury najwyższe i Jak widać 
są one rządu 900 K, Dowodzi to słuszności sformułowanego w 
pierwszym rozdziale wniosku o nieprzydatności modeli 
hydrodynamicznych do opisu zjawiska odwrotnej kumulacji. 
Temperatury maksymalne rządu 900 K, a tym bardziej 
temperatury średnie rządu ^ 400 -r 600 K są bowiem znacznie 
mniejsze od temperatury topnienia miedzi, która wynosi 1356 K. 
Vyłątek może stanowić jedynie wstępna faza napędzania wkładki. 



nazwana fazą hydrodynamiczną, gdyż mimo niskich temperatur jest 
wtedy spełniona relacjai ^ 

Nie mieszczący sią również w ramach hydrodynamicznej 
teorii efekt wyrównywania prądkości "zbitka" i "strumienia" 
ilustrują od strony ilościowej zmiany w czasie prądkości 

i . JTa wykresie tym naniesiono również ustaloną 

wartość prądkośct pocisku zmierzoną eksperymentalnie. Jak widać 
rozbieżności wyników teoretycznych 1 eksperymentalnych w 
zakresie tego parametru są rządu 7%. Dokładność tego rządu- Jest 
zupełnie wystarczająca z punktu widzenia aplikacji praktycznych 
omawianej metody, Z.drugiej strony wobec złożoności problemu 1 
umowności szeregu modelowych założeń nie mają Już głąbszego 
uzasadnienia próby uzyskania wląkszej dokładności. 

Rozkłady gąstoścl średniej wzdłuż przekrojów C-C i D-D' 
sygnalizują ciekawe i dość złożone zjawisko rozwarstwiania sią 
płyt metalowych napądzanych ciśnieniem produktów detonacji. 
Zjawisko to obserwuje sią dla płyt. o dostatecznej grubości, 
które rozdzielają sią na dwie lub nawet wlącej cząści wzdłuż 
płaszczyzn równoległych do ich powierzchni. Efekt ten, znany 
Już dawno od strony eksperymentalnej £2,141 Jest w dalszym 
ciągu przedmiotem badań teoretycznych, w tym również 
wykorzystujących metody fizyki komputerowej <np,prace 
£46,1131). 

V ramach niniejszej pracy zagadnienie rozwarstwiania sią 
wkładek kumulacyjnych w trakcie ich napądzania było Jednym z 
podprobłemów, który musiał zostać rozwiązany w taki sposób, aby 
mimo efektów zniszczenia materiału w obszarach rozwarstwień, 
móc kontynuować obliczenia dla dalszych faz zjawiska kumulacji. 
Rozwiązanie takie umożliwił oczywiście omówiony wcześniej model 
powstawania i wzrostu objątości szczelin. ¥ modelu takim obszar 
rozwarstwienia bądzie obszarem bardzo silnego, lokalnego 
wzrostu objątości szczelin (spadku gąstości średniej) z 
praktycznie zerowymi składowymi tensora naprążśń < J 2® * 
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Jak jut wspomniano efekty takie sygnalizują jut rozkłady 
gestc-śct wzdłui przekrojów C-Ć' .1 X>-D* ( g <0.7 Jest w 

pr .-ybliienlu równoważne >, ale dla ich zrozumienia 

trzeba prześledzić dokładniej czasowo ■" przestrzenny rozwój 
tego zjawiska.tf tym celu można posłużyć sie ilustracjami 
zamieszczonymi na rys.25. 

Na pierwszej ilustracji C i =8.75 ps) widoczny jest 
wydłużony obszar zniszczonego materiału usytuowany praktycznie 
równolegle do powierzchni wkładki. Jest to wi^c typowy efekt 
rozwarstwiania-si^ wkładki, V tym miejscu mogą zrodzicS sie dwa 
pytania, a mianowicie : dlaczego taki efekt w ogóle powstaje 
craz dlaczego powstaje w obszarze zaznaczonym na rys.25 ? 

V pracy [143, w oparciu o przybliżenie akustyczne 
pokazano, że rozwarstwienie występuje wtedy, gdy Impuls 
ciśnienia generowany na kontakcie; produkty detonacji - wkładka 
1 propagujący się we wkładce ma charakterystyczny kształt 
nny. Otóż, na odległości równej dwóm grubościom 
wkładki musi on zmniejszać się o wielkość przekraczającą 
progowa wartość na zniszczenie materiału . Wówczas w 

odległości X od powierzchni swobodnej wkładki naprężenia 
ro*.1 a^ 1 na czole tali odciążenia mogą przekroczyć wartość 6^ 

i rozpocznie sie proces rozwarstwiania. Zakładając, że Impuls 
ciśnienia postać: 

P-P»{''-4J <4.4) 


mc^na zgodnie z pracą (143 podać następując^ oszacowanie 
wielkości X; 




JL <=' 
z f>® 


(4.5) 


Napre.<,enie ^ można byłoby utożsamiać z wprowadzonym w modelu 
tworzenia się szczelin naprężeniem progowym , ale okazało 

się z obliczeń, że dynamicznie naprężenie to może być nawet 
kilkakrotnie przekroczone ^szczególnie na początku procesu gdy 
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Ilyo. 25 . Efelct rozwarstwiania się wkładlti łctmiulacyjnoJ w procesie Jej napędzania. 










Vc - v^< > Dlatego też dla oceny wleli:o?^ci X lepiej przyjąć 
zgodnie z prac^ [46;. że dla miedzi ^ 3x10^^ dyn/cm®. 

Aby ccenld charakterystyczna dla naszego przykładu wielkości 
1 to naniesiono i\& rys. 35, kształt Impulsu ciśnienia pA= PaW 
zdjętego w punkcie A* który zaznaczono strzałką na pierwszej 
ilustracji. Przyjmując szacunkowo na podstawie tego wykresu 
po ^ 3x10'’ dyn/cro® ł ^ 2,5x10''^ s dostajemy z <4.5): 

X ^ 5x10-® cm. 

Oszacowar.le to dość dobrze zgadza sią z wynikami symulacji 
komputerowej, gdyż Jak to widać na rys.25, w dość znacznym 
otoczeniu punktu A rozwarstwienie zaczyna sie na odległości 
X Cs 6xi0*® cm. Nie można jednak z pomocą wzoru <4.5> 
wytłumaczyć wszystkich obserwowanych efektów. Ma on bowiem 
bardzo przybliżony charakter 1 nie uwzględnia na przykład 
faktu, że samo przekroczenie naprężenia 6*^ nie oznacza jeszcze 
zniszczenia ośrodka. Do tego potrzebne Jest Jeszcze odpowiednio^ 
długie utrzymanie naprężeii przekraczających wartości progowe, ; 
aby objętość szczelin mogła wzrosnąć do poziomu V<ł >(*Vc)|^^. 

Ton efekt oraz złożony obraz falowy w otoczeniu podstawy 
wkładki są cdpowiedzialne za charakterystyczy kształt 
zniszczonych obszarów w tym rejonie. Analogiczny obraz 
znisz<:zon uzyskano również na bazie symulacji komputerowej w 
pracy [1143 1 potwierdzono go eksperymentalnie. Natomiast w 
ci::.r.tralnych obszarach wkładki rozwarstwienie Jest natychmiast 
ii widowano przez efekt pogrubiania sie wkładki (symetria 
cylindryczna) w wyniku ruchu masy w kierunku osi <U<0>. 

Na dalszych ilustracjach ( t = 16 ps 1 t = 30 ps) 
cami eśzczonych na rys. 25 pokazano ewolucję w czasie kształtów 
zniszczonego obszaru. Charakterystyczne Jest przy tym 
pcdawienle daifczych znl^^zczeń w wyniku radialnego 

rozciągania powierzchniowych warstw wkładki, które nałożyło się 
na wcześniejszo zniszczeni© w obszarze rozwarstwienia- Kształt 
zniszczonych obszarów pokazanych na tych Ilustracjach tłumaczy... 
r-.t w pełni rozkłady gęstości pokazane na rys. 23. 


Omawiając podstawowe równania problemu zasygnalizowano, że 
opis efektów lepkich w dynamicznie deformowanych metalach Jest 
w dalszym ciągu sprawą otwartą 1 to nie tylko dlatego, że 
stosuje się różne modele teoretyczne, np,Maxwella, Binghama, 
model dyslokacyjny i inne. Okazuje się bowiem, że nawet w 
ramach Jednego modelu oceny, ilościowe efektów lepkich podawane 
przez różnych autorów mogą się różnić nawet w zakresie kilku 
rzędów wielkości. Przykładem może tu być najprostszy loodel 
Maxwella ze współczynnikiem lepkości ||^=const, V literaturze 
można znaleźć przykładowo następujące oszacowania tego 
współczynnika: 

10 P w pracy C 433, 2x10* P w C115,1163, 3x10* P w t1173 oraz 
10^ -j- 10® P w pracach t40, 118, 1193, 

V tej sytuacji przyjęto również w pierwszym etapie pracy 
stały współczynnik lepkości, dobierając Jego wartości z zakresu 
10* -s* 10® P. Okazało .się, że dla dużych -i- 10® P, 

modelowane teoretycznie pociski miały wyraźnie mniejszą długość 
1 większą średnicę niż w eksperymentach. V miarę zmniejszania 
10* ” 10* P relacje te poprawiły się, ale równocześnie 
zaczął się pojawiać efekt oderwania ramion wkładki od Jej 
centralnej części, czego nie obserwowano w eksperymentach. 

Szczegółowa analiza tych wyników wskazywała wyraźnie, że 
lepkość ma istotny wpływ na parametry wybuchowo kształtowanych 
pocisków, ale założenie =const jest zbyt dużym 
uproszczeniem. Dlatego też w dalszych badaniach uwzględniano 
zależność współczynnika lepkości od temperatury, która 
teoretycznie powinna być zależnością najsilniejszą, 2 powodu 
braku w literaturze ilościowych ocen dotyczących tej zależności 
przyjęto wstępnie, że funkcja F(T} Jest liniowa (patrz wzór 
2,32>, a stąd: 

r -i. - (.1- «•' 

Parametr Tm został określony przy omawianiu modelu Steinberga. 
Optymalną wartość współczynnika dla miedzi ustalono na 


















poziomie 10^ P. Uzyskano przy tym znaczni© bliższe w stosunku 
do eksperymentu wyniki, co dowddzlio poprawności sformułowanych 
wniosków 1 przyj<ątego kierunku badań. 

Aktualnie najlepsze wyniki osiągnięto wprowadzając w 
miejsce zależności liniowej <4,6) z^ależność wykładniczą: 

exp(-4^) [-1+ bp(-||)*/^ - h (T- 500)] <4.7) 

Postać funkcji przyjęta analogicznie Jak dla stopionej 

miedzi [120}, a w ramach badań własnych ustalona tylko 
współczynnik =5 x 10“® P, Trzeba Jednak powiedzieć, że 

różnice wyników dla ^ issiją Już Jakościowego 

charakteru, ale dotyczą raczej szczegółów związanych z 
kształtami badanych pocisków. Bóżnice te można ocenić porówując 
wyniki przedstawione na rys.26 i 27, gdzie w Jednakowych 
skalach przestrzennych 1 Jednakowych odstępach czasu 
przedstawiono wyniki dla współczynników lepkości określonych 
wzorami <4.6) 1 (4.7). 

Jednocześnie na rys,28 przedstawiono fotografię wybuchowo 
ukształtowanego pocisku, który był przedmiotem badań 
teoretycznych,. Dla porównania zamieszczono obok najlepszy 
aktualnie wynik teoretyczny w tej samej skali przestrzennej 
(ostatni kadr z rys.27), 

Z rys,28 widać, że zgodność wymiarów 1 kształtów pocisków Jest 
zupełnie wystarczająca z punktu widzenia zastosowań 
praktycznych omawianej metody, - z 

V kontekście tych rozważań można zwrócić Jeszcze uwagę na 
dość nieoczekiwane zastosowanie zjawiska "odwrotnej kumulacji”. 
Otóż, można Je wykorzystać, tak Jak zrobiono to w niniejszej 
pracy, do określenia niemierzalnych bezpośrenlo w eksperymencie 
charakterystyk materiałowych, Jak np. . Kod 

komputerowy w połączeniu z odpowiednimi rejestracjami 
eksperymentalnymi (fotografiami w świetle widzialnym lub 
Jeszcze lepiej rentgenowskim) staje się praktycznie nowym 
narzędziem w dziedzinie badah własności materiałów. 
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Bys.26.EwolucJ 
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Rys. 28 Pardvmanie vr/niJcóv7 eksperyjusTitalnych i 

teoretycznych dotyczące ostatecznych kształtów 
wybuchowo uformowanego pocisku. 

4,2. Proces tworzenia klasycznego strumienia kumulacyjnego. 

Bazując na opraćowanym modelu matematyczno-fizycznym oraz 
dof^wladczeiiiach. natury numerycznej zdobytych, przy badaniu 
zjawiska ‘‘odwrotnej kumulacji" mo±na było przystąpić do 
zweryiikowania przydatnoścl metody w badaniach innych 
z jawisk. Qozvwtocie » problematyka pracy wyiBagałat aby w 
pierwszej kolejności podjąć próbę zamodslowanla procesu 
tworzenia Gią klasycznego struiiiienia kumulacyjnego. 

podobnie jak poprzednio przedmiotem badań był układ o 
skali laboratoryjnej, średnica wkładki wynosiła 3 cm, grubość 
0.1 cm, a połówkowy kąt rozwarcia wkładki oC przyjęto równym 
30^. Rbv/riież, podobnie jak w poprzednim przykładzie zastosowano 
do badań termoplastyczny materiał wybuchowy. 

Kształt badanego ładunku oraz szczegółowy 'w:.;gląd 
pokrywającej go sieci numerycznej pokazano na rys,29. 












z kolei na rys.30 zilustrowano proces detonacji materiału 
wybuchowego (punkt inicjacji detonacji zaznaczono strzałką). 
Widoczne są na nim również początki procesu deformacji wkładki 
i generacji strumienia kumulacyjnego. 

ITa ostatnim kadrze zamieszczonym na rys, 30 pokazano sytuacją, w 
której wpływ ciśnienia produktów detonacji na wkładką jest już 
pomijalnie mały i dlatego wyeliminowano je z dalszych obliczeh. 

Deformacją wkładki oraz proces tworzenia slą strumienia 
kumulacyjnego zilustrowano na rys.31 1 32.. 

Oczywiście z rysunków tych można określić tylko czasowo - 
przestrzenne zmiaąy kształtów badanego układu. Ze wzglądu na 
to, że skala tych rysunków nie pozwala prześledzić szczegółowo 
procesu deformacji sieci numerycznej, zamieszczono na rys.33 
przykład powiąkszonej sieci ilustrujący charakter Jej 
deformacji. 

Po zapoznaniu sią z czasowo-przestrzennyM 
charakterystykami geometrycznymi zjawiska przedstawionymi na 
rys. 29 -5- 33 można przejść dp omówienia jego podstawowych 
charakterystyk fizycznych. Najciekawsze charakterystyki dotyczą 
w tym przypadku strumienia kumulacyjnego i Jego parametrów 
takich, jak; prądkość, temperatura, gąstość i masa. Rozkłady 
tych parametrów wzdłuż ośl "z** zilustrowano na rys. 34 w fazie, 
w której strumień jest już całkowicie ukształtowany i ulega 
Jedynie dalszemu wydłużaniu slą skutkiem istnienia gradientu 
prądkośći poosiowej, 

Z analizy zamieszczonych na rys.34 wykresów wynikają 
nastąpujące wnioski: 

a/ Osiowa składowa prądkości masowej ma praktycznie 

liniowy rozkład wzdłuż osi *'z‘\ a prądkość zbitka jest 
wyrównana i ma typową wartość ^^0.5 km/s. 
b/ Temperatura czoła strumienia wskazuje, że tworzy je 
stopiona miedź (T > 1400 K>. V zasadniczej cząścl 
strumienia temperatura Jest Jednak niższa od temperatury 
topnienia 1 zgodna z eksperymentalnie określonymi 
wartościami “ 1200 1300 K [ 2J . 









c/ Rozkład gęstości wzd2u± strumienia wskazuje, że utworzony 
Jest on z litego metalu < 0), z wyjątkiem ewentualnie 

niewielkiego obszaru na samym czole strumienia. W fazie , 
dla której sporządzono te wykresy, nastąpiło już wyraźne 
oderwanie sie strumienia od zbitka, Świadczy o tym silny 
spadek gęstości średniej w obszarze oddzielającym strumień 
od zbitka ( Vo iCf ). Zjawisko to jest 

udokumentowane.eksperymentalnie E21. 
d/ łlasa strumienia Mj = H/(z) jest liczona od czoła 

strumienia do danego punktu na osi "z*\ Masa aktywnej 
cząścl strumienia <z SJ 10 -t- 21 cm> stanowi około 25 % 
całkowitej masy wkładki Mc . 
e/ Parametry takie, Jak: prędkość zbitka,prędkość czoła 

strumienia, temperatura średnia strumienia, gęstość i masa 
aktywnej części strumienia są typowymi wartościami dla 
klasycznych ładunków kumulacyjnych omówionych w pracy [2]. 
świadczy to już o ogólnej poprawności otrzymanego 
rozwiązania, które później będzie Jeszcze dodatkowo 
skonfrontowane z wynikami własnych eksperymentów. 

Ra rys.31 l 32 widoczny jest bardzo wyraźnie wzrost w 
czasie kąta rozwarcia wkładki. Począwszy od trzeciego kadru na 
rys,31 średnia wartość tego kąta Jest juź większa od 90 To 
ciekav/e zjawisko Jest wynikiem bardzo słabego obciążenia 
podstawy wkładki Impulsem ciśnienia generowanym w tym obszarze. 
Riewtelka ilość materiału wybuchowego t boczne fale odciążenia 
powodują bardzo szybki zanik tego Impulsu w wyniku czego 
podstawa wkładki uzyskuje bardzo małą prędkość, zarówno 
radialną jak i osiową. Rie bez znaczenia jest również fakt, że 
początkowy kąt rozwarcia wkładki o^=30«>, jest kątem dość dużym 
dla tego typu eksperymentów, W układach praktycznych zjawisko 
to w takiej skali z reguły nie występuje, gdyż naddatek 
materiału wybuchowego i osłony mogą znacznie zwiększyć impuls 
ciśnienia w tym obszarze. 

V omawianym przykładzie nie wystąpił tak wyraźnie jak 
poprzednio efekt rozwarstwiania się Wkładki, Jest ona bowiem 
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2.5 rara cieii-ssa, a ponadto prawie d.wuł:rotnie maleje w tym 
przypadku amplituda odbitej lall detonacyj tiej Cspadek kata 
rozwarcia wkładki cL Ju± z prsyblijioaeęa oszacowania C4.5> 
wynika, ie mamy w tej sytuacji X *?• L, gdzie L oznacza grubość 
wkładki, a rozwarstwienie mażą wystąpić dla X<L. Warunek X<L 
rcoźe być spełniony tylko w pewnym otoczeniu podstawy wkładki, 
gdyż spada tam również czas trwania impulsu . Z rozwiązania 

komputerowego wynika, że istotnie, w otoczeniu podstawy wkładki 
CT i-3 ^ nastąpiło przekroczenie granicznego 

naprężenia ó"^ i rozpoczął proces zniszczenia materiału. 

Jednakże niewielka grubość wkładki powoduje dodatkowo skrócenie 
czasu, dla którego zachodzi ffl > iw tym konkiretnym 

przykładzie czas ten był zbyt krótki, aby spełnił się warunek 
zniszczenia >(Vc)j^^- Maksymalna wartość objętości szczelin 

Vc nie przekroczyła w żadnym punkcie wkładki wartości około 
fVc)Kr . Kożna więc rozważania te podsumować stwierdzeniem, 
źG w omawianym przykładzie nie nastąpiło rozwarstwienie się 
wkładki, ale jedynie pe//ne osłabienie materiału w otoczeniu Jej 
podstawy. 

Bazując na rozważanym X)rzykład2ie symulacji komputerowej 
procesu tv/orzenia się strumienia kumulacyjnego można dokonać 
ciekav;Gj l bardzo konkretnej weryfikacji wyników otrzymywanych 
w oparciu o przedstawioną w rozdziale I teorię hydrodynamiczną. 

Oczywiście, podstavrawym warunkiem umożllwlającyia jej 
zastosowanie, jest znajomość parametrów napędzonej wkładki, 
tzn. rozkładów prędkości masowej i kąta jaki tworzy ona z osią 
symetrii. V naszym przypadku Informację te najprościej uzyskać 
v; oparciu o rozv?iązania komputerowe w raoiLencie, w którym 
zakończył się Już proces napędzania wkładki. V takim właśnie 
moji^^ncle <t-13 na rys. 35 przedstawiono kształt wkładki wraz 

z naniesionym rozkładem prędkości masowej 

(W wybranych punktach wzdłuż linii 1>, V oparciu o te 
inlorniacje można określić średni kąt oC jaki wkładka tworzy z 
osią oraz przybliżony rozkład prędkości W wzdłuż 
linii l. 
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Przyjęto zatem, zgodnie z oznaczeniami występującymi we wzorze 
<1.16>! VVt>-3 km/s, ł1d=0,86 oraz o4 =55'=^* 

Dane te vłystarczają do zbudowania dla *^>13 ps rozwiązania w 
oparciu o teorię hydrodynamiczną, które można będzie 
skonfrontowaó z rozwiązaniem numerycznym. Wyniki obu rozwiązań 
przedstawiono n- rys,36, Widzimy, że zarówno kształty strumieni 
jak i rozkłady prędkości masowej pozostają ze sobą w bardzo 
dobrej zgodności. 

Oznacza to, ża w badaniach procesów tworzenia się strumieni 
kumulacyjnych podstawowe znaczenie ma znajomość parametrów 
napędzonej wkładki .kumulacyjnej, Jeśli bowiem nie dysponujemy... 
możliwością otrzymania pełnych rozwiązań problemu to możemy 
ograniczyć się tylko do jego fazy wstępnej, w trakcie której 
odbywa się napędzanie wkładki. Dalsze rozwiązanie można już 
konstruować w oparciu o teorię hydrodynamiczną. Podejście takie 
zastosowano w pracy [261, gdzie do określenia stanu napędzonej 
wkładki wykorzystano kod typu Lągrange*a, korzystając z faktu, 
że w fazie- napędzania wkładką! jej deformacje są stosunkowo 
nlowieikie. 

Należy jednak podkreślić, ż© uzyskanie tak dobrej 
zgodności wyników <rys.36> było możliwe dlatego, że badany 
przykład w czasach 13 i: ^37 ps spełnia warunki 

stosowalności przybliżenia hydrodynamicznego <obszar 1 na 
rys.6), Ze wzoru <1.13) wynika bowiem, że próces jest cały czas 
poddżwiękowy - M 0,5, a efekty wytrzymałościowe określone 
warunkiem <1,40> zaczną być istotne dopiero dla elementów 
wkładki, które w chwili t-13 pś miały promień większy od około 
l cm. Ponieważ elementy ta zaczną tworzyć strumień dla 
t >37 p-s, to stąd wynika właśnie czas obowiązywania 
przybliżenia hydrodynamicznego *- t<37 ps. Wykresy zamieszczane 
na rys,36 dotyczą zatem maksymalnych czasów, dla których 
przybliżenie to można stosować. Rozwiązanie numeryczne 
przedstawiono w chwili wcześniejszej o 4 ^s, tak aby czoła obu 
strumieni miały to samo położenie na osi z. Przesunięci© to 
wynika z faktu, że rozwiązanie hydrodynamiczne rozpoczyna się w 


































Błoniencie gdy w rozwiązaniu numeryczny aa czoło strumienia 
Crys.35> jest już ukształtowane i odległe od punktu S o ponad 
2 cm. Jest to wynikiem róvmoczesności procesów napądzania 
wkładki 1 tv^orzenia si^ czoła strumienia* Mieuwzglądnienie tego 
faktu w przybliżeniu hydrodynamicznym nie prowadzi do dużych 
błędów, a wynika to stąd, źe efekty te dotyczą pomijalnie 
małych mas tworzących czoło strumienia. 

Porównanie otrzymanych rezultatów z dostępnymi 
Informacjami literaturowymi oraz teorią łhydradynamiczną daje 
Już określoną podstawą do stwierdzenia iż są one poprawne 
zarówno ad strony jakościowej Jak i ilościowej. Tym niemniej, w 
celu uzyskania bezpośredniego dowodu tej poprawności wykonano 
odpowiednią serią eksperymentów. 

V ich wyniku otrzymano serią fotografii strumienia w wybranych 
chwilach czasu. Kilka tych fotografii zamieszczono na rys.37, 

Sa rys.37 a/ widoczne Jest Jedynie spalające sią w atmosferze 
czoło strumienia. Fotografią tą można wykorzystać w charakterze 
bazy do wyznaczania prądkości czoła strumienia. Na nastąpnej 
fotografii r- 37 b/ widoczny jest już sam strumień kumulacyjny 
oraz jego silnie świecące czoło. ¥ oparciu o te dwie fotografie 
wyznaczono prędkość czoła strumienia, która podobnie jak w 
rozwiązaniu .lUJim^^cznym, wynosi około 6 fcm/s. Ną bazie 
fotografdii 3T b^/ naożna również w przybliżeniu ocenl,ć średnicą 
struml'siviLćs kŁimniiIL 2 GcyJ nego w określonej odległości od Jego czoła. 
Rozbieżności między teorią a eksperymentem w tym zakresie są 
mniejsze niż wynosi błąd odczytu średnicy strumienia na 
zamieszczonej fotografii. Fotografia 37 c/ pokazuje fazą silnie 
wydłużonego strumienia i początki jego fragmentacji, a 37 d/ 
fazą zaawansowanej fragmentacji. Do problemu tego powrócimy 
jeszcze w tym rozdziale. 

Prace eksperymentalne potwierdziły więc bezpośrednio (w 
zakresie mierzonych parametrów) poprawność uzyskanych rozwiązań 
1 tym stwierdzeniem można już zakończyć on^wlanie wyników 
uzyskanych dla modelowego układu służącego do wytwarzania 
klasycznych strumieni kumulacyjnych. 
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Rys.37,Fotografie strumienia kumulacyjnego wytworzonego w 
badanym układzie w kolejnych chwilach czasuJ 
a/ -17.5 psj b/ -30,5 psj 
c/ -48.5 psj d/ -66.5 ps. 
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Pozytywne wyniki tsstów i weryflkac.ji eksperymentalnych 
dla dotychczas ainGv/lcnych układGw kumulacyjnych stały sIę 
bcdzcein da podjęcia prób zaJtucdolowania układów produkowanych w 
skali przemysłowej. Wbrew pozorom przejście od układów o 
charakterze " laboratoryj nym** do układów *'przemysłowych” nie m 
tylko formalnego charakteru zmiany wymiarów tych układów. 

W układach tych należy sie bowiem liczyć z koniecznością 
zaiiiadeiQvrania szeregu dodatkowych efektów, które nie 
występowały w omówionym układzie eksperymentalnym. 
Scharakteryzujmy krótko najważniejsze z tych efektów. 

1. y dotychczasowych badaniach przyjmowano, ż© detonacja była 
inicjowana na osi symetrii układu, a co za tym idzie front 
fali detonacyjnej miał kształt sferyczny. Tymczasem w 
praktyce często wykorzystuje się odpowiednie przesłony, za 
pomocą których można sterować propagacją frontów fal 
detonacyjnych tak, aby materiał wybuchowy bezpośrednio 
oddziałytnający z wkładką był inicjowany na pewnym 
promieniu wokół osi. V modelu teoretycznym fronty tak 
inicjowanych fal detonacyjnych będą więc miały kształt 
toroldalny, a nie sferyczny, 

2. Materiał wybuchowy Jest z reguły umieszczany w 
odpowiednich osłonach, które są Jego mechanicznym 
zabezpieczeniem. Jednakże osłony te oddzlaływując z falami 
detonacyjnymi zmieniają czasowo przestrzenny kształt 
impulsu ciśnienia obciążającego wkładkę i muszą być 
uwzględniane w obliczeniach. 

3. Układy ”przemysłowe” charakteryzuje bardzo duży stosunek 
średnicy wkładki do Jej grubości, który może osiągać nawet 
wartości rzędu '*^10®, Stwarza to poważne problemy natury 
technicznej związane z budową odpowiedniej dla takich 
przypadków sieci numerycznej. Ponadto często stosuje się w 
praktyce wkładki o zmiennej grubości, co dodatkowo 
komplikuje te problemy. 

4. Omówione w punkcie i nieosiowe inicjowania detonacji ma na 
celu zwiększenie impulsu ciśnienia działającego na 
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wkładkę. Dotyczy to głównie wierzchołkowych części 
wkładki, z których powstaje czoło strumienia. V związku z 
tym należy się liczyć z możliwością pojawienia się 
przepływów naddżwlękowych w fazie tworzenia się czoła 
strumienia. Jak wiadomo z rozważań przeprowadzonych w 
rozdziale I, przepływom takim powinna towarzyszyć 
fragmentacja strumieni. Kod komputerowy musi więc być 
przygotowany na eliminowanie z obliczeń takich fragmentów, 
które przestały oddziaływać fizycznie z pozostałą częścią 
strumienia. 

Efekty omówione* w punktach 1,3 1 4 dają się zamodelować 
stosunkowo prosto, gdyż wymagają jedynie odpowiedniej adaptacji 
omówionych wcześniej elementów metody numerycznej. Natomiast 
modelowanie wpływu osłon na przebieg zjawiska kumulacji jest 
problemem nowym 1 dlatego trzeba go omówić nieco dokładniej. 

V ścisłym sformułowaniu należałoby modelować dynamikę 
ruchu osłon analogicznie jak zrobiono to dla wkładki 
kumulacyjnej. SkomplIkowałoby*to jednak znacznie kod numeryczny 
i wymagałoby dalszego zwiększania czasów obliczeń oraz obszarów 
zajmowanej pamięci. Problemów tych można uniknąć modelując ruch 
osłony w sposób przybliżony, tym bardziej, że można ją wkrótce 
po napędzeniu wkładki syyellminować z dalszych obliczeń, gdy jej 
wpływ na modyfikację impulsu ciśnienia staje się pomijalny. 

V przykładzie, który będzie zilustrowany na rys.38,39 1 40 
osłonę zamodelowano ciągiem mas skupionych, których ruch 
opisano następującym układem równań; 
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B ilość elementów na jakie podzielono osłonę, 

- masa L-tego elementu (stała), 

- składowe prędkości elementu odpowiednio wzdłuż osi 
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^ - aktualny promień i długość elementu, 

f*-- kąt Jaki twarzy element z osią z, 
pj_ - ciśnienie produktów detonacji na L-ty element. 
Warunek brzegowy na kontakcie produkty detonacj1 - osłona 
modeluje się analogicznie Jak warunek na powierzchni kontaktu 
produkty detonacji - wkładka kumulacyjna. Każdy punkt produktów 
detonacji oddzlaływujący z osłoną znajduje najbliższe sobie 
masy skupione reprezentujące tę osłonę. W oparciu o ich 
prędkości można zaraodelować warunek swobodnego poślizgu tak 
Jak omówiono to w poprzednim rozdziale, natomiast dla każdego 
elementu (masy skupianej) reprezentującego osłonę, w oparciu o 
aktualne ciśnienia najbliższych mu punktów reprezentujących 
produkty detonacji, wylicza się średnie ciśnienie 
wymuszające ruch elementów - wzory C4.8> -5- <4.9). Aby uniknąć 
przypadkov^ch fluktuacji prędkości elementów wskazane jest 
wprowadzenie pewnej korekty uśredniającej i wygładzającej 
rozkłady prędkości elementów sąsiednich. Może ona mieć 
przykładowo następującą postać: 
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Współczynnik liczbowy A dobiera się eksperymentalnie 
(typowo A,'5fe0.01-^0.03>, 

Na rys.38,39 i 40 pokazano przykład rozwiązania 
numerycznego dla klasycznego ładunku kumulacyjnego, w którym 
v/ystępują omówione wyżej elementy charakterystyczne dla 
ładunków stosowanych w praktyce. Promień, na którym inicjowana 
jest detonacja zaznaczono strzałkami na rys,38. Elementy osłony 
ilustrują znaczniki złażone z układów trzech punktów. Na rys.40 
przedstawiono sytuację, w której można już wyeliminować z 
dalszych obliczeń produkty detonacji wraz z elementami osłony. 
Nie zamieszczano już ilustracji rozwiązań dla późniejszych faz 
zjawiska, bo różnią się one od układu modelowego tylko skalą 
czasowo - przestrzenną. 
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Rys.38, 39, 40.Deformacje sieci numerycznej pokrywające 

materiał wybuchowy, wkładkę kumulacyjną i osłonę, t 






































Rys.41.Rozkłady prędkości masowej , średniej gęstości g i 
temperatury T wzdłuż osi symetrii tworzonego w 
przepływie naddźwiekowym czoła strumlerxla. 


Z punktu widzenia fizyki zachodzących procesów ciekawsza 
jest faza początkowa, w której tworzy slą czoło strumienia, Vr 
fazie tej mamy 'bowtara do czynienia z omówionym w rozdziale I 
przepływem naddżwl^kowym (skutkiem silnego obciążenia 
wierzchołka wkładki przy nieoslowej inicjacji detonacji)* 
Maksymalną liczby Macha w omawianym przykładzie można oszacować 
na około 1.7. Na ryś,41 przedstawiono kształt czoła strumienia 
oraz rozkłady wzdłuż osi symetrii; prędkości masowej , 

średniej gęstości g i temperatury *T . Widać, że czoło 
strumienia tworzy ośrodek, który można traktować tak jakby , 
składał sl^ z praktycznie niezależnych odłamków - Vc >CVt)Kr' 
Ze zmian gęstości średniej wynika, że proces fragmentacjt 
pogłębia Sie w miarę zbliżania sią do granicy czoła strumienia* 

Zgodnie z tym co powiedziano w rozdziale I prędkość 
2c ^ 8 km/s (dla miedzi) powinna oddzielać obszar 
przepłyynj naddżwiekowego od poddż wiekowe go, a wlec strumienia 
rpzfragmentowanego od litego. Jak widać z rys.41 oba warunki: 

i 2c wskazują praktycznie na ten sam punkt 

strumienia, Świadczy to zarówno o poprawności przeprowadzonych 
w rozdziale I ocen jak i dobrej ocenie krytycznej objętości 
szczelin. . 

Kończąc rozważania dotyczące klasycznych ładunków 
kumulacyjnych należy zwrócić Jeszcze uwagę na proces 
przemieszczania sle 1 wydłużania strumienia pod wpływem 
gradientu prędkości poosiowej. Procesowi temu towarzyszą bowiem 
narastające niestabilności przewężeniowe, które po określonym 
czasie prowadzą do rozpadu strumienia na szereg podążających za 
sobą niezależnych elementów. Początek zjawiska można 
zaobserwować na fotografii 37c/, a na fotografii 37d/ oddzielne 
elementy są już widoczne bardzo wyraźnie. 

Niestety proęes ten nie może być należycie zamodslowany w 
oparciu o tworzącą strumień sieć punktów obliczeniowych, gdyż w 
fazie fragmentacji nie Jest już ona wystarczająco gęsta. Aby 
móc śledzić rodzące sie niestabilności trzeba byłoby zage^*^^^^ 

Ją przynajmniej o rząd wielkości (długość elementów, na które 















rozpada si^ strujmień Jest rz^du 3 do 5 Jego średnic). Badanie 
tego zjawiska mogłoby być realne tecłmlcznle gdyby po 
zakończenia procesu tworzenia sie struinieaia pokryć go nową, 
interpolowaną siecią punktów obliczeniowych, kosztem 
wyeliminowanego z dalszych obliczeń zbitka. 

Podjecie w przyszłości tego typu badań wydaje sie 
nieuniknione, gdyż? 

1. Prace teoretyczne dotyczące tego problemu są bardzo 
nieliczne. Rezultaty wynikające z teorii małych zaburzeń 
podana w pracy C41] - dla przybliżenia hydrodynamicznego 
oraz w 1121] - dla ośrodka sprężysto -- 
plastycznego.Natomiast w pracy 1122) omówiono próbą 
komputerowego modelowania procesu rozwoju niestabilności 
przeważeniowych, tworząc zagęszczoną sieć pokrywającą 
pewien wycinek strumienia. 

2. Wynik tego typu badań ma decydujące znaczenie przy ocenie 
głębokośól przebicia konkretnych układów kumulacyjnych 
[2,22]. Nie wnikając szczegółowo w te problemy można 
powiedzieć, że zjawisko fragmentacj1 strumienia Jest 
odpowiedzialne za Istnienie optymalnej odległości wkładki 
od pancerza, zapewniającej maksymalną głębokość przebicia 
l14). Gdyby zjawiska tego nie było, głębokość przebicia 
musiałaby ciągle wzrastać w miarę oddalania wkładki od 
pancerza. 

4.3.Przebicie osłony ciałem odkształcałnym. 

Globalnym sposobem oceny efektywności ładunków 
kumulacyjnych są parametry drążonych przez nie otworów w 
różnego rodzaju osłonach (pancerzach). Parametry te porównuje 
się dla różnych konstrukcji odnosząc Je z reguły do średnicy 
ładunku kumulacyjnego. Przy przebijaniu osłon strumieniem 
kumulacyjnym takim podstawowym parametrem Jest głębokość 
przebicia, a przy wykorzystaniu wybuchowo kształtowanych 
pocisków również średnica tworzonego krateru. Możliwość 


teoretycznej oceny tych parametirów wymaga opanowania 
problematyki symulacji komputerowej procesów przebijania osłon 
ciałami odkształcalnymi. Problematyka ta nie stanowi Jednak 
przedmiotu niniejszej pracy i nie będzie szerzej rozwijana. V 
rozdziale tym zostanie Jedynie zademonstrowane, źe opracowana 
metoda rozwiązywania problemów kumulacyjnych, może być z 
powodzeniem zastosowana również do badania efektów przebijania 
osłon. Stwarza to możliwość zbudowania w przyszłości kodów, 
które obejmowałyby Jeszcze szerszy niż obecnie kompleks 
zagadnień od zainicjowania detonacji począwszy, a na przebiciu 
osłony skończywszy. 

Na rys.42 pokazana przykład modelowania procesu przebicia 
miedzianej płytki o grubości 1 cm, miedzianym walcem o średnicy 
2 cm 1 długości 2 cm, napędzonym do prędkości 2 km/s. 

Wybór miedzi na materiał zderzających się ciał wynikał 
oczywiście z faktu dobrej znajomości charakterystyk 
wytrzymałościowych dla tego metalu. Natomiast kształty l 
wymiary tych dal przyjęto kierując się wyłącznie względami 
natury technicznej, gdyż opracowana metoda pozwala modelować 
zderzenia ciał o różnorodnych kształtach. 

Na dwóch ostatnich kadrach rys.42 widoczna Jest 
f ragmentacj a ” powłoki** tworzącej się na granicach tarczy i 
uderzającego ciała, Zasady na Jakich wykonywano te operacje są 
oczywiście takie same Jak w poprzednich przykładach. Różnica 
polega tylko na tym, źe w omawianym przykładzie, po 
wyeliminowaniu części **powłokl'*, dalsze obliczenia prowadzone 
są równolegle dla dwóch nleoddzlaływujących z sobą obszarów - 
pozostałości po tarczy 1 uderzającym walcu. Zabieg taki nie był 
poprzednio wykonywany. 

Porównanie pokazanych na rys.42 deformacji zderzających 
się ciał z przykładowymi rozwiązaniami literaturowymi (np.CGl), 
[70}> wskazuje na niewątpliwą poprawność uzyskanego 
rozwiązania, przynajmniej od strony jakościowej. Głębsze 
analizy dotyczące zgodności ilościowej uzyskanych rozwiązań z 
teoretyczno - eksperymentalnymi doniesieniami literaturowymi 









Rys.42.Deformacje sieci numerycznej ilustrujące proces 

przebijania miedzianej tarczy, miedzianym walcem o 
prędkości 2 km/s. Odstęp czasowy miedzy kadrami-13.6 ps. 
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nie bądą przeprowadzane, gdyż przykład ten ma tylko 
demonstracyjny charakter. Tym niemniej wydaje siQ celovłym 
sformułowanie kilku uwag dotyczących tego problemu. 

Z wykresów zamieszczanych w pracy i 123] można wnioskować, 
że średnica krateru w tym przykładzie powinna być około dwa 
razy większa od średnicy uderzającego walca. Tego. rsqdu 
średnica krateru Jest widoczna na 3 i 4 kadrze rys.42. Ua 
dalszych kadrach widać jednak powolne zwiększanie <się tej 
średnicy. Zasadniczą przyczyną tego zjawiska Jest to, że tarcza 
ma stosunkowo małą średnicę i w wyniku przebicia pierścieh Jaki 
z niej pozostał zaczął w całości ekspandować radialnie. 

Oprócz parametru Jakim jest średnica krateru 
najwygodniejsze do interpretacji 1 porównah są parametry fali 
uderzeniowej i struktura obszaru powstającego w centrum 
zderzających sią ciał. Ifa rys. 43 przedstawiono rozkłady 
gęstości ^ i prędkości masowej wzdłuż osi układu w kilku 

wybranych chwilach czasu. 

Dla czasu "t = 0.91 ps z zamieszczonych W 5 ^kresów można 
odczytać parametry fali uderzeniowej rozprzestrzeniającej się w 
obu kierunkach osi ”z** od oznaczonej literą 3 powierzchni 
kontaktu uderzającego walca z tarczą (przyjęto, że czas liczony 
Jest od momentu zderzenia się ciał>. Ponieważ uderzający walec 
miał prędkość 2 >m/s to jest rzeczą oczywistą, że prędkość 
masowa fali uderzeniowej musi wynosić 1 km/s i taka też 
prędkość ustala się w otoczeniu powierzchni iS . 2 adiabaty 
uderzeniowej miedzi C2] można odczytać, że gęstość za frantem 
takiej fali powinna wynosić 10.9 g/czo® co również bardzo dobrze 
zgadza się z wartością zaznaczoną na wykresie. Parametry za 
frontami fal uderzeniowych są więc określane bardzo precyzyjnie. 
Pewne zdziwienie może natomiast budzić duża szerokość 
numerycznie modelowanego frontu fali uderzeniowej. Z rys.42 
widać, że tarczę podzielono zaledwie na 7 "komórek” wzdłuż osi, 
a typowe rozmycie frontu przez pseudolepkość wynosi od 3 do 5 
komórek. Tłumaczy to tak duże rozmycie frontu fali uderzeniowej 
i jednocześnie Jest cenną charakterystyką metody 
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nu me rycz Be J , która mimo \^^rowadzenla 'ZQvmetrzBeJ dyfuzji 
numerycznej nie rozmywa frontów falowych w stopniu większym niż 
inne metody wykorzystujące pseudolepkośó. Oczywiście w 
obliczeniach o charakterze użytkowym należałoby w takim 
przypadku znacznie zagęścić sieci przestrzenne^ aby zmniejszyć 
strefy rozmycia frontów falowych 1 tym samym zwiększyć 
dokładność metody obliczeniowej. 

tf chwili czasu tT - 3,7 jis w tarczy propaguje się Już fala 
odciążenia, a swobodna, zewnętrzna powierzchnia tarczy rusza z 
prędkością około 1.3 km/s. Vynik ten na pierwszy rzut oka 
wydaje się być wyraźnie zaniżony. Vtadomo bowiem, że w tym 
zakresie parametrów swobodne powierzchnie metali uzyskują 
prędkości bliskie podwojonej prędkość! masowej fali 
uderzeniowej C21, a więc w naszym przypadku powinna to być 
prędkość około 2 km/s. Okazało się jednakże, że poosiowa 
prędkość masowa około 1 km/s tylko w pewnym otoczeniu 

powierzchni zderzenia S • Fala uderzeniowa wygenerowana w 
otoczeniu tej powierzchni jest bowiem falą ekscentryczną i 
ąuaslsferyczną, a więc jej parametry będą się obniżać w miarę 
oddalania się frontu od powierzchni S . W otoczeniu 
powierzchni swobodnej prędkość masowa na czole fali 
uderzeniowej spada skutkiem tego efektu do około 0,7 km/s i 
dlatego w fali odciążenia uzyskuj© ona prędkość zaledwie 
1,3 km/s. Analogicznie tylna, swobodna powierzchnia walca nie 
zostaje wyhamowana do prędkości 1 km/s ale Jak widać z rys.43 
do około 1,5 km/s, Dla ścisłości należy dodać że mogą 

być obarczone pewnym błędem wynikającym ze zbyt dużej strefy 
rozmycia frontu fali uderzeniowej. 

2 wykresów dla późniejszych chwil czasu widać, że od 
strony obu powierzchni swobodnych nastąpiło rozwarstwieni© 
ośrodka. Wystąpiło ono jednak bardzo blisko tych powierzchni, 

CG przy niedostatecznej gęstości sieci numerycznej objawiło się 
jedynie powstaniem obszarów □ niskiej średniej .gęstości, bez 
możliwości dokładnego rozpoznania Ich struktury. V miarę upływu 


161 











czasu widać również sxc3pniovry zanik prcce-saw falowych i 
związane z t y Jii wyr ć v/ny wa nie s i ^ r oz ki adó w pr ą dko -ś o 1 laa s owe j , 
Kończąc rozważania zawarte w tym rozdziale nałoży 
podkro-.!ić, źe zagadnienie przebijania osłon Jest w raiiiacli 
mechaniki ośrodków .ciągłych jednym z na.jdłużej i najsilniej 
rozwijanych kierunków wykorzystujących metody symulacji 
komputerowych. V dalszym ciągu ukazuje si^ wisie prac z tej 
dziedziny^ np. [ 45, 1244-129, 131] » przy czym aktualnie zwi-aca sią 
w nich uwagę głównie na doskonalenie 1 weryfikację opisów 
własności materiałowych. 

Omówione przykładowe rozwiązanie dowodzi, ze opracowana i 
opisana w pracy metoda symulacji komputerowej również nadaje 
się do badań tego typu problemóv/. Rozszerza to z Jednej strony 
zakres metod komputerowych stosowanych w tej dziedzinie, a z 
drugiej strony umożliwia podjęcie w przyszłości kompleksowych 
badań problemów kumulacji i przebijania osłon, w ramach Jednego 
matematyczno - fizyczno - numerycznego modelu. 


Rozdział 5 


yKIOSKl 

V pracy zebrano i przedstawiono materiał dotyczący 
opracowanego i zweryfikowanego w praktyce sposobu komputerowego 
modelowania zjawi-sk kumulacyjnych. Na tle krótkiej 
charakterystyki tych zjawisk oraz aktualnego stanu wiedzy na 
ten temat omówiono komleksowo proces modelowania, poczynając od 
modelu matematyczno - fizycznego, poprzez budowę odpowtedniego 
kodu numerycznego, a kończąc na przykładach konkretnych 
rozwiązań. ¥ rozdziale ninie.1szym zebrane zostaną najważniejsze 
vmioskl wynikające z przedstawianej pracy. Mają one charakter 
pewnego podsumowania osiągniętych wyników t jednoczesnego 
wskazania działań, których podjęcie umożliwiłoby w przyszłości 
dalszy, krajowy rozwój tych badań o bardzo dużym znaczeniu 
praktycznym. 

1. Głównym celem prowadzonych aktualnie eksperymentów w 
dziedzinie kumulacji Jest budowa układów o bardzo wysokich 
parametrach i dużej skuteczności (głębokość przebicia, 
objętość krateru, itp,), Zbudowanie układów mogących 
sprostać współczesnym wymaganiom wyłącznie na drodze' 
eksperymentalnej staje się stopniowo nierealne, zarówno ze 
względu na konieczną w takim przypadku, dużą liczbę 
eksperymentów Jak t kłopoty z Ich interpretacją. Pracom 
takim potrzebne Jest więc -silne i kompleksowo opisujące 
eksperyment wsparcie teoretyczne, 

2. Na obecnym etapie prac eksperymentalnych takiego wsparcia 
nie może Już dać im klasyczna teoria kumulacji. 2akres JeJ 
stosowalność! powoduje bowiem, że nie znajduje ona 
zastosowania w badaniach procesów wybuchowego 
kształtowania i napędzania pocisków, a przy analizie 
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zjawisk twar:3'3nla si^ strumieni kujnułaoyjnych JsJ 
zastosowanie wymaga uprzedniej znajomości parametrów 
napędzonej wkładki. Jest to oddzielny i znacznie 
trudniejszy problem niż samo zastosowanie prostych wzorów 
klasycznej teorii. 

3. Od połowy lat 70-tych poJav/laJą sIę informacje 
literaturowe o próbach homputeroweso modelowania procesów 
iruaiulacyj nych. V badaniach te^o typu można uwzglądnlaó 
nlestacjonarnośó i wielowymiarowość przestrzenna zjawisk 
oraz realne warunki brzegowe 1 własności materiałów. 
Kierunek ten otwiera wiec nowe, szerokie laożliwości 
kompleksowego badania zjawisk fizycznych w całej ich 
złożoności. 

4. Informacje literaturowe dotyczące problemów modelowania 
komputerowego zjawisk kumulacyjnych mają Jednak charakter' 
przykładowych 1 fragmentarycznych wyników dotyczących 
wybranych wariantów układów kumulacyjnych. Ule ma 
możliwości uogólnienia tych .yynlków, a wlec i skorzystania 
z nich w .,elu wsparcia własnych prac eksperymentalnych. 
Dlatego też. aby wykorzystać wszystkie możliwości Jakie 
stwarza modelowanie komputerowe, konieczne Jest opanowanie 
określonych metod 1 posiadanie własnych kodów 
komputerowych, które umożliwiają na bieżąco dostąp do 
wszystkich interesujących nas informacji o badanym' 

Układzie, 

V pracy pr.,,edstawlono własny, w znacznym stopniu 
oryginalny,matematyczno - fizyczno - numeryczny model 
opisujący procesy detona.jJ i i dynamicznej deformacji ciał 
stałych, 26 szczególnym uwzglądnlenlem zjawisk 
kumulacyjnych. Wykorzystując zjawiska "klasycznej" 1 
■odwz-otnej" kumulacji zweryfikowano ten model od strony 
Jakościowej i tlośclov,eJ. Aktualnie oslągnląta zgodność 
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wyników teoretycznych i eksperymentalnych jest zupełnie 
wystarczajćica z punktu widzenia prac doświadczalnych. 
Ogólnie inożna Ją ocenić na dokładność rzędu kilku procent 
w zakresie wszystkich mierzonych parametrów. 

Przedstawiona metoda daje szerokie możliwości modelowania 
złożonych zjawisk z dziedziny mechaniki ośrodków ciągłych, 
a w szczególności: 

zjawisk opisywanych równaniami gazodynamlki, 
hydrodynamiki, teorii sprężysto - plastyczności lub 
sprężysto/lepko - plastyczności, 

procesów nlestacjonarych i dwuwymiarowych przestrzennie, 
warunki brzegowe mogą być stawiane na krzywoliniowych 1 
ruchomych powierzchniach, 

rozwiązywane mogą być równocześnie sprzężone problemy 
opisywana równaniami gazodynamlki i mechaniki ciał 
stałych, 

możliwe JfcSt modelowanie bardzo dużych deformacji ośrodka, 
łącznie z jego rozwarstwianiem się lub fragmentowaniem na 
niezależne części. 

V ramach dotychczas prowadzonych prac, omówiona metoda 
symulacji komputerowej znalazła konkretne zastosowanie 
przy badaniach następujących zjawisk fizycznych: 
wybuchowe kształtowanie i napędzanie pocisków, 
tworzenie się strumienia kumulacyjnego, 
przebijanie osłon ciałami odkształcalnymt. 

Pokazano przy tyra, że metodę tę można z powodzeniem 
stosować do ładunków o charakterze przemysłowym 
<uwzględniając dodatkowo: realne relacje wymiarów, a 
zwłaszcza średnicy do grubości wkładki, osłony, nieosiowe 
inicjowanie detonacji, itp.). 

Analiza dotychczas uzyskanych rezultatów wskazuje', że 
sposród zjawisk ogólnie scharakteryzowanych w punkcie 6, 
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V-' plerv5roj kolejności można by toby podjąć badania 
r.3stepujących problemów: 

analiza śtabiinoóci i ira^mentacj1 strumieni 
ku cj 13 c y j Ti yc )i , 

P^^^bijanie osłon ciałami o nieregularnych kształbach 
<np. w:/bychDwo formowane pociski lub strumienie 
kumulacyjne), a w dalszej ^jperspektywie połączenie 
zaradnieh kumulacji 1 przebijania osłon, 
generacja supersilnych fal uderzeniowych w układach 
stoikowych lub walcowych, 

analiza- procesów napędzania i stabilnoiści koncentrycznych 
otoczek <linerów>, 
wybuchowe łączenie metali» 

wybuchowe prasowanie proszków w układach wielowarstwowych. 

charakter metody numerycznej stwarza szerokie 
iłoźllwoScl dalszego jej udoskonalania i badania Jej 
właściwości. Zwłaszcza, że jest Ona źródłem szczególnie 
wielu problemów natury ideowej i technicznej, które można 
rozwiązywać na wiele sposobów. Najważniejsze r nich to; 
wybór punktów sąsiednich, 
eliminacja różnych typów niestabilności, 
dobór zewnętrznej dyfuzji numerycznej, 

metody obliczeh efektów lepko — plastycznego płjmi^cia i 
tworzenia sl^ szczelin, Itd. 

Ar punktu widzenia teorii metod różnicowych bardzo ciekawe 
mogłyby być pogłębione badania warunków stabilności 
omówionej metody, gdyż analiza przeprowadzona w pracy 
, n.i,ej1 miała tylko Jakościowy charakter. 

również podkreślić, że opracowana metoda numeryczna 
ca/kowicie otwarta na modyfikacje matematyczno - 
znogo opisu zjawisk, w tym również na dalszą jego 
Iowę. w zależności od potrzeb, kierunku dalszych prac lub 
■wt.i nowych Informacji literaturowych. 



10, Róvmalegle z badaniem określonych zjawisk fizycznych 
prowadzone były prace zmierzające do określenia własności 
niateriałov/ych badanych wkładek. V szczególności dotyczyły 
one ustalenia v/s pół czynników i zależności funkcyjnych 
opisujących efekty lepko/plastycznego płynięcia 1 
tv/orzeala sle szczelin wraz z określeniem ich wpł 3 rwu na 
własności wytrzymałościowe. 

Tego typu badania własne były niezbędne, gdyż doniesienia 
literaturowe na ten temat są nieliczne i znacznie 
zróżnicowane u różnych autorów, Należy liczyć sie z tyra, 
że w dalszym, ciągu bedzie pojawiać sl^ konieczność ich 
kontynuowania, zwłaszcza jeśli bedzie stą rozszerzał 
zakres prowadzonych prac. 

11, Opracowana metoda numaryczna, podobnie Jak metoda PIC. 
uraożliv/iać powinna wzglądnle proste modelowanie deformacji 
ośrodków wielofazowych. Sprawdzenie 1 potwierdzenie tych 
możliwości stworzyłoby szansy na podJe'^ie badań 
przebijania pancerzy o konstrukcji warstwowej Lub 
generc^cii strumieni kumulacyjnych z wykorzystaniem wkładek 
blme taiic z nyo h. 

12, Kodęlovr.anie komputerowe zjav/isk fizycznych może być 
skutecznym narzędziem wspierającym prace eksper 3 ?Ji!entaln.e. 
Jednoczę.śnie jest ono w znacznym stopniu od tych pr.BC 
zależne t to z dwóch co najmniej powodów. Wyniki 

eksDerymentalne są bowiem podstawą jakościowej i 
ilościowej weryflk,acjl uzyskiwanych rozwiązań oraz dają 
możliwość v/^/znaczenia charakterystyk materiałowych 
niedostępnych w literaturze, 

V badaniach zjawisk kumulacyjnych szczególnie przydatne są 
dwa rodzaje diagnostyk; fotografia w świetle widzialnym t 
fotografia rentgenowska. Dotychczasowe prace teoretyczne 
były konf rontcv/ane tylko z fotograr ‘ sirai wykonanymi w 
świetle widzialnym. Dalszy postęp tych prac będzie zaleśny 











w duiryni sToprau od rozwoju diagnostyki rentgenowskiej » 
niogacel dostarczyć bardziej szczegółowych, inlarmacji o 
badanych zjawiskach. Diagnostyka ta jest stosowana i 
doskonalona w liczących sie w :twieęle ośrodkach badawczych 
za;’mu lacych si-ą zbliżoną problematyką [1321, 

K3',omj,asi do badań własności dynamicznie deformowanych 
materiałów coraz czó^ciej wykorzystuje si^ metody 
1nterferometryczne. Umożliwiają one wykonanie bardzo 
precyzyjnych pomiarów prędkości powierzchni swobodnych 
deiorroowanych 1 napędzanych ciał. Uzyskane w ten sposób 
w;^/r.lkl są konfrontowane z wynikami symulacji komputerowej, 
uiuożl 1 wia ląc wyznaczenie poszukiwanych charalcterystyk 
materiałowych. Eksperymenty wykonuje sią przy tym w taki 
sposób, aby przy ich interpretacji mogły być 
yr/korzystywane Jednowymiarowe kody numeryczne. Daje to 
iTiOżilwość pi'zeprowadzerila dużej liczby testów 
komputerowych w stosunkowo krótkim czasie, a wląc czyni ‘ 
proces określania własności materiałów szybkim 1 

t ywnym. yetoda interferometryczna <wraz z odpowiednim 
t s t o w V m k G.ii e m n i oe r y c z n y m) powinna r ówn i e ź Jak 

Vv .:ej w:.:bo,pacie zestaw diagnostyk stosowanych w 
b a d a n .1 a c r* kra ! □ y/’/c h . 
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534.222.2.001.5 Mechanika. Detonacja . Badanie . Modelowanie . 

539.3 Mechanika ciał stałych. Odkształcenia. 

519.66B _ Matematyka komputerowa . Algorytmy. _ Pol. _ 

Jach K.: Modele nie komputerowe zjawisk kumulacyjnych 
Warszawa: WAT 1990, 179s.rys.tab.bibliogr. poz.132 

W pracy przedstawiono kompleksowy fizyczno-matematyczno-numeryczny model służący do opisu 
dużych.dynamicznych deformacji metali, detonacji materiałów wybuchowych oraz oddziaływania 
produktów detonacji na ciała stałe. Podstawowymi zjawiskami do opisu których zastosowano ten 
model była tzw.klasyczna lub "odwrotna" kumulacja oraz dodatkowo efekt przebijania osłony 
ciałem odkształcalnym. Praca zawiera kompleksowy opis oraz przykładowe rozwiązania powyższych 
problemów (z wykorzystaniem aktualnej wiedzy o własnościach dynamicznie deformowanych materia¬ 
łów i własnych badań teoretyczno-eksperymentalnych). Rozwiązania teoretyczne skonfrontowano 
z wynikami eksperymentów wykazując ich zgodność w stopniu całkowicie wystarczającym dla prak¬ 
tycznych zastosowań modelu. Przeprowadzono pomyślnie zakończone próby bezpośredniej adaptacji 
modelu do optymalizacji ładunków kumulacyjnych o charakterze przemysłowym. 

Model fizyczno-matematyczny oparty jest na równaniach teorii lepkoplastyczności w sformułowa¬ 
niu Hohenemsera-Pragera i uzupełniony o półempiryczne równania stanu oraz zależności opisujące 
zmiany dynamicznej granicy plastycznego płynięcia.modułu ścinania i współczynnika lepkości w 
funkcji temperatury, ciśnienia, gęstości i deformacji plastycznej (model Steinberga). Ponadto 
w obszarach niszczenia struktury metali korzystano z fenomenologicznego modelu opisującego pow¬ 
stawanie i wzrost objętości szczelin oraz ich wpływ na charakterystyki wytrzymałościowe. 
Przykładowe rozwiązania prezentowane w pracy oraz dane liczbowe wynikające z własnych badań 
dotyczą wkładek kumulacyjnych wykonanych z miedzi - podstawowego do chwili obecnej materiału 
konstrukcyjnego w tej dziedzinie. 

Procesy detonacji opisano klasycznymi równaniami gazodynamiki oraz półempirycznymi równaniami 
stanu. Front fali detonacyjnej aproksymowano powierzchnią silnej nieciągłości. 

Od strony matematycznej model jest zespołem kilkunastu nieliniowych równań różniczkowych cząst¬ 
kowych w trzech zmiennych niezależnych (czas i dwa wymiary przestrzenne) uzupełnionych związka¬ 
mi algebraicznymi opisującymi własności materiałów oraz odpowiednimi warunkami początkowo^brze- 
gowymi. 

Dla rozwiązania tak sformułowanego problemu zaproponowano, przebadano i opisano w pracy ory^T*^ 
nalny algorytm symulacji komputerowej typu "punktów swobodnych". Umożliwia on rozwiązywanie 
szerokiej klasy problemów z zakresu mechaniki ośrodków ciągłych, gdyż pozwala na: prowadzeniś 
obliczeń w warunkach skrajnie dużych deformacji (łącznie z fragmentacją), stawianie warunków 
brzegowych na ruchomych i krzywoliniowych powierzchniach, zszywanie rozwiązań na kontakcie 
metal-ciało stałe, korzystanie z różnych modeli opisujących własności ciał stałych, cieczy, 
gazów, produktów detonacji itd. 
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